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ABSTRACT 



Ocean acoustic tomography permits the mapping of various properties of a body 
of water through indirect means. The technique utilizes travel time variations for an 
acoustic signal to determine the structure of the ocean medium via inverse mathematical 
methods. The scale of any tomography experiment is fundamentally limited by the 
signal to noise ratio at the receiver. Through the use of a near vertical acoustic array, 
normal mode modeling of the local environment and a modal beamformer, array gains 
are possible which greatly extend the maximum separation between source and receiver. 
Additionally, the technique provides temporal resolution of the modal components of 
the arriving signal. 

A time domain modal beamformer for a near vertical acoustic array has been 
developed. It has realized a noiiinal array gain of 6 dB for the Heard Island 

Experiment vertical array deployed off California. The primary obstacle to the 
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INTRODUCTION 



I . 



A. THESIS SUMMARY 

The objective of this thesis is to develop a software 
package to beamform acoustic signals used in ocean acoustic 
tomography. The goal of tomography signal processing is the 
precise measurement of acoustic travel time. Sound speed is 
a well understood function of temperature, pressure and 
salinity. Therefore, a fluctuation in acoustic travel time is 
indicative of changes in the environment through which the 
acoustic energy has passed. Once the travel time fluctuations 
are known, inverse mathematical techniques can be used to 
infer various properties of the ocean medium. [Ref. 1] 

One source of travel time measurements is through the use 
of explosive devices. Although such signals are easy to 
generate, they may not provide the required precision owing to 
dispersion of the various frequency components in the 
impulsive signal. Additionally, such signals are difficult or 
impossible to replicate. [Ref. 2] 

A better method that has been employed in recent years is 
the use of maximal -length sequences . The technique utilizes 
a pseudorandom, binary sequence as the phase modulation for an 
electronically generated bandpass signal. Maximal-length 
sequences are well suited to travel time measurements by 
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virtue of their deterministic nature, autocorrelation 
properties and simplicity. 

The goal of the work described in the following pages is 
to provide sufficient gain to permit the recovery and decoding 
necessary for accurate travel time measurements for acoustic 
paths of extreme length (in this case 10,000 nautical miles). 
Specifically, the programming package should be able to: 



• Exploit the spatial structure of the incident wave field 
to discriminate among signals arriving from different 
directions . 

• Utilize information regarding time varying array tilt and 
depth to estimate the array geometry. 

• Provide a stable virtual array by using the array geometry 
and modal structure of the immediate environment. 

• Distinguish among the various modal components of the 
target signal. 

• Provide both time domain and frequency domain analysis. 

• Remain sufficiently flexible so as to permit the 
processing of arrays of differing construction in future 
tomography experiments. 

• Provide some measure of fault tolerance with respect to 
the receiving array. 



The remainder of this thesis is structured as follows: 



• Chapter II. 

• Chapter III. 

• Chapter IV. 

• Chapter V. 



Acoustic Wave Propagation Theory 
Modal Beamforming 
The Heard Island Array 
Results and Conclusions 
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Chapter II presents an introduction to acoustic wave 
propagation. The approach taken is that of normal mode 
propagation in a range independent channel. 

The third chapter is an introduction to the concept of 
modal beamforming. It reviews frequency domain and time 
domain beamf ormers . The algorithm utilized in this study is 
derived. 

Chapter IV describes the construction and subsequent 
deployment of the receiving array used in this experiment. 
Data acquisition and preprocessing of the acoustic signal are 
discussed . 

The final chapter presents a description of the array 
dynamics encountered and frequency domain output from the 
software developed. Additionally, proposals for system 
integration and performance improvements are detailed. 

Appendix A contains the source code which has served as 
the body of this work. In addition to the beamformer, various 
utility programs are included which provide for array geometry 
data reduction. 

B. THE HEARD ISLAND EXPERIMENT 

The performance of the software package is evaluated on a 
raw data set (vice synthetic data). The data was acquired 
during the Heard Island Experiment which occurred during the 
period January 23 - February 2, 1991. The purpose of the 
experiment was to determine the reliability of global acoustic 
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paths for tomographic analysis. Specifically, it is desired 
to establish the viability of these transmission paths for a 
proposed multi-national attempt to detect a decreasing trend 
in acoustic travel time. Such a change would indicate an 
overall warming along the path, providing the first convincing 
evidence to the existence of global warming [Ref. 1]. 

The signals were transmitted from the vicinity of Heard 
Island in the southern Indian Ocean. The site was selected 
because it is central to a web of open water paths extending 
through all the worlds oceans as shown in Figure 1.1. 




Figure 1.1: Heard Island Raypaths 



The transmitting ship was the R/V Cory Chouest. The 
transmitter utilized a ten element vertical array with a 
maximum of five elements active at any time. A nominal source 
level of 209 dB re 1 pPa at 1 meter was realized. Individual 
projectors employed were HLF4LL very low frequency sources. 
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A variety of signal types were sent, including continuous 
wave and maximal -length sequences. A carrier frequency of 57 
Hz was chosen both for its low absorption losses and the 
ability to distinguish it from the 50 and 60 Hz frequencies 
generated by power plants world wide. 

The Monterey component of the Heard Island Experiment 
involved a collaboration between the Naval Postgraduate 
School, the Monterey Bay Aquarium Research Institute, the 
Massachusetts Institute of Technology, Woods Hole 
Oceanographic Institution, SAIC and the Moss Landing Marine 
Laboratory. The R/V Point Sur deployed a 32 element vertical 
array approximately 70 nautical miles south-west of Monterey 
Bay . 
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II. ACOUSTIC WAVE PROPAGATION THEORY 



A. ACOUSTIC PROPAGATION IN AN OCEAN WAVEGUIDE 
1. The Inhomogeneous Wave Equation 

The homogeneous linear wave equation [Ref. 3], 



V*p = 



dt^ 



( 2 . 1 ) 



incorporates two implicit assumptions. The first restricts 
the sound speed to a constant value in the region of interest. 
Following the development of Coppens [Ref. 4], if the sound 
speed is now permitted to vary in space (but not time), then 
c=c(x,y,z) may be substituted into the wave equation without 
loss of generality. 

The second implicit assumption restricts the use of 
the wave equation to sound fields which are free of sources. 
Before including a source term, the selection of a specific 
coordinate system is appropriate. 

Given that sound in the ocean does not spread 
spherically in a free field, but radially with upper and lower 
boundaries, implies the use of a cylindrical coordinate 
system. To further simplify the problem to one sufficient for 
this work, the sound field shall be assumed to exhibit radial 
symmetry about the source, and sound speed variations shall be 
restricted to depth, c(z). Based on these restrictions. 
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application of the Lapacian yields the homogeneous, range 
independent wave equation, 






r dr\ dr c^(z) 

The inclusion of a source term requires that the 
inhomogeneous wave equation reduce to homogeneous form in 
regions which are free of sources of acoustic energy. This 
requirement implies that the delta function for a point source 
located at r=0 and z=z^ have the form [Ref. 3] 






( 2 . 2 ) 






6 (r) 6 (z- z^) , 



(2.3) 



where r is the spherical radius vector, r is the cylindrical 
radial coordinate and z is the cylindrical depth coordinate. 

Inclusion of this terra in the wave equation yields the 
Helmholtz equation for a monofrequency point source of unity 
amplitude. 



IJ.1 






r dr\ 


i drj 


a^2 \ciz)j 



--|6 (r) 6 (z-z^) 



(2.4) 

2 . Solutions to the Inhomogeneous Wave Equation 

Having adopted a cylindrical coordinate system, 
equation (2.4) can be solved by separation of variables. 
Furthermore, the complete solution can be treated as a linear 
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combination of normal modes, 



p(r,z,t) = i?„(r) Z^(z) . ( 2 . 5 ) 

m 

The normal modes (Z„) form an orthonormal set of 
eigenfunctions in z which satisfy the sourceless Helmholtz 
equation, 



dz^ 



(0" 



V (z) 




0 , 



( 2 . 6 ) 



subject to the appropriate surface and bottom boundary 
conditions, and the orthonormal condition. 



fZniz)Zjz)dz = b^, (2.7) 

where is the eigenvalue for the m^^ eigenfunction (or normal 
mode) and is the Kronecker delta function . 

A closed form expression for the R„(r) may be obtained 
by substitution of equation (2.6) into the inhomogeneous wave 
equation (2.4), multiplication of all terms by integrating 
over z and application of the normalization condition (2.7). 
The result, after manipulation, is 



\ _d 
I dr 



dr 






-j6(r)Zjz,), 



( 2 . 8 ) 



the inhomogeneous Bessel's equation. It has the known 
solution 



K(- 






(2) 



(K.r) 



(2.9) 
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where is the Hankel function of the second kind and 
order zero. Substitution of this form into equation (2.5) 
yields 

pir.z.t) = Z„(z) Z„(z„) _ (2.10) 

m 

Closed form solutions for the Z„(z) are either 
difficult or impossible for all but select cases. 
Fortunately, efficient numerical algorithms exist which 
provide for the rapid solution of the depth dependent 
functions. [Ref. 5] 
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III. MODAL BEAMFORNING 



A. PRELIMINARY CONCEPTS 

Propagating waves can be modeled as functions of both 
space and time. Consequently, the attributes of such models 
can be used to extract information from real wave fields. 
Beamforming exploits the temporal and spatial characteristics 
of a specified environment to enhance a particular aspect of 
the wave field while attenuating undesirable components. Such 
an operation can be loosely termed constructive reenforcement 
of the desired signal or noise rejection. The information 
exploited in this development includes the modal structure of 
the immediate environment and the instrument data supplied by 
the array itself. 

B. TIME DOMAIN BEAMFORMING 

The acoustic signal received at a given hydrophone can be 
expressed as 

p{t,x,y,z) = S(t,x,y, z) +N(t,x,y, z) , (3*1) 

where p is the pressure at the hydrophone, S is the desired 
signal and N is the local noise field. All are functions of 
time and hydrophone location. 

If one now considers the signal received at an array of 
hydrophones, equation (3.1) takes the form 
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Pit) = Sit) *Nit) , 



(3.2) 



where p, S and N are all N x 1 vectors representative of the 
individual array elements. The location of individual 
elements is implied by the vector index for a given 
hydrophone . 

Classical (or plane wave) beamforming incorporates the use 
of a complex steering vector (E) to impose a phase shift on 
the individual elements in order to enhance the sensitivity of 
the array to signals propagating in a specific direction, 



E = 



J^2 



(3.3) 



e 






The output signal from such an array is 

bfpjt) = iSit) ^Nit)) . (3.4) 

The subscript (pw) indicates that the beamformer assumes the 
presence of plane waves. Additionally, the implicit 
assumption is made in the above expression that the 
autocorrelation function of the noise field will be 
identically zero for any two (distinct) hydrophones. The 

success or failure of the beamformer will be a direct 
reflection of the validity of these assumptions. 
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If the desired signal is not monof requency, but of finite 
bandwidth, simple phase shifting may not be suitable. This is 
particularly true if the target signal is carrying information 
such that the communication bit time is comparable to, or 
greater than the time taken for the wave to propagate across 
the array. Under such circumstances, the beamformer would 
impose phase shifts in distinct communication bit segments 
simultaneously. The resulting distortion to the bandpass 
signal is unacceptable in most communications applications. 
Beamforming under these conditions requires the application of 
true time delays vice phase shifts. 

Assuming an array of spatial extent such that time domain 
beamforming is required implies a summation of the form 

N 

bfp^it) = ^ (5( t-tn) + AT( t-T„) ) , (3«5) 

J3 = l 

where n is the index used to describe the array elements and 
represents the steering delay applied to the n^^ element 
[Ref. 6]. This arrangement places no additional restrictions 
on the target signal characteristics or array construction in 
a given physical system. 

C. THE MODAL BEAMFORMER 

The signal model employed above is incomplete, in that it 
does not account for the nonuniform distribution of acoustic 
energy in the sound channel. Vertical arrays inherently 
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sample the modal structure of the immediate environment. 
Therefore, one needs to take the information provided by 
normal mode theory into account when the acoustic array spans 
a region of nonuniform energy density. 

The appropriate refinement to the complex steering vector 
which incorporates the modal nature of the signal arrival is 



Z2(zi) 

2 ^ (Zj) Z 2 (Zj) 



Z2(z„) 



Z^(z2)e^*- 



Z^(z^) 



(3.6) 



where the m^^ column of the matrix is the generalized steering 
vector for the m^^ normal mode. The individual steering 
vectors now include amplitude factors (Z^(z^)) that reflect the 
value of the eigenfunction at the various hydrophone depths 
[Ref. 7]. The n^^ row contains the values of those steering 
vectors at the depth of the n^^ hydrophone. The phase terms 
(?nm) account for phase shifts required for the beam steering 
applicable to all modes (as required by array geometry) and 
phase delays among the various modes due to their individual 
propagation speeds. Specifically, the take the form 



% ' < 3 - 7 ) 

where 0^ is determined by the hydrophone location and is 
determined by the differences in propagation speeds over the 
transmission path for the normal modes present. 
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Modifying the steering matrix to implement time delays 
(vice phase shifts) is straight forward. A beamformer could 
be implemented by summing the output of individual hydrophones 
(p(t)) at this point. However, to do so implies a 
deterministic evaluation of the a„. Such an evaluation would 
require near perfect knowledge of the conditions over the 
transmission path prior to commencement of the experiment. 
Given that the aim of tomography is the mapping of these 
conditions, such an approach is inappropriate. 

An alternate approach is to utilize the columns of the 
matrix individually to discriminate among the modal components 
of the signal. Since all components of a given modal steering 
vector share a common value of these phase terms may be 
discarded from the individual elements of the steering vector. 
As a result of this simplification, the modal steering vector 
is completely determined by the (obtained from normal mode 
modeling of the immediate environment) and the desired look 
direction (determined by array geometry). Having established 
the methodology, and replacing the phase weights with steering 
delays, the output of the modal beamformer is 

73 = 1 

where the subscript (m) indicates that the desired output is 
the component of signal representing the m^^ normal mode. 
Expanding this form to explicitly include the terms 
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representing the signal and noise components of the acoustic 
field yields 

( 3 . 9 ) 

n=l n=l 

where and represent the amplitudes of signal and noise 
at the n^^ hydrophone. The delay and weighting of the noise 
received at the hydrophones has no undesirable effects. If 
the noise between any two hydrophones was uncorrelated prior 
to the operation, then it will remain so after weighting. 
Again, the degree to which the ambient noise field is 
uncorrelated will ultimately reflect upon the performance of 
the beamformer. 
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IV. THE HEARD ISLAND VERTICAL ARRAY 



A. ARRAY CONSTRUCTION 

The array deployed for this experiment consisted of 32 
hydrophones deployed vertically, each having a nominal 
sensitivity of -170 dB re 1 V/pPa. The element spacing was 45 
meters, with hydrophone number one occupying a design depth of 
345 meters. The nominal design depth for hydrophone number 32 
was 1740 meters. 

Instrumentation on the array consisted of two sensors. 
The upper sensor was located 4.0 meters above hydrophone 
number one. It recorded ambient pressure, temperature, tilt 
and current velocity. The lower sensor was located 5.0 meters 
below hydrophone number 20. It recorded tilt (including 
direction), pressure, temperature and conductivity. This 
sensor appears to have suffered a casualty during array 
deployment which rendered its tilt data suspect. 

Flotation was provided by one main syntactic float and 28 
syntactic football floats. The design called for the main 
float to reside at a nominal depth of 230 meters. 
Approximately half of the football floats were submerged. 
This arrangement rendered the array neutrally buoyant, thus 
providing a measure of isolation from surface wave effects. 
The array was directly tethered to the R/V Point Sur with 
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floatation devices on the surface portion. Each float was 
equipped with a flasher, a flag and a radar reflector. The 
float nearest the ship had the addition of a radio beacon. 
The procedure called for the ship to remain approximately 1200 
meters from the position immediately above the array. During 
the experiment the R/V Point Sur was configured to remain 
quiet and dead in the water, except as required to keep clear 
of the array. Figure (3.1) illustrates the general 
configuration . 




Ballast ^ 



Figure 3.1: Heard Island Receiving Array 
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B. INSTRUMENT DATA 



Ideally the receiving array, when deployed, is oriented 
vertically. However, due to currents, the requirement to 
maintain control of the research vessel and other factors 
beyond the control of the research team, the array is often 
tilted from the vertical. Given the length of the array, a 
1.0° tilt results in a 24 meter horizontal displacement 
between hydrophones one and 32. This displacement is nearly 
one (carrier frequency) wavelength. Additionally, the 
direction in which the array tilts affects the component of 
horizontal displacement which is collinear to the propagation 
direction. Horizontal displacements perpendicular to the 
direction of propagation are acceptable if one assumes that 
the wave field is fairly uniform is this direction. 
Therefore, the array steering utilized in this work seeks only 
to correct for the component of displacement in the direction 
of signal propagation. 

The speed with which the wave field propagates is required 
prior to the calculation of steering delays. The appropriate 
speed to use in this case is the group speed of the carrier 
frequency for the mode in question, defined as 



c 



ffin 




( 3 . 10 ) 



where Cg^j^ is the local group speed for the m^^ normal mode, o> 
is the radian frequency and is the appropriate eigenvalue. 
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The amplitude weight applied to a given hydrophone is a 
function of the mode number sought and hydrophone depth. The 
mode shapes (Z^(z)) are required prior to beamforming. The 
numerical algorithm used to obtain the eigenvalues and 
eigenfunctions for this work was developed by Chiu and Ehret 
[Ref. 5]. It utilizes finite differencing of the propagation 
vectors for a given frequency and sound speed profile to 
calculate the modal structure of the local environment. The 
sound speed profiles were measured during the conduct of the 
experiment by the research team aboard the R/V Point Sur. 
Figure 3.2 is the sound speed profile used for this analysis. 
Figures 3.3 through 3.5 are the output of the normal mode 
model. Amplitude weights for the array elements were 
calculated from these data sets. 
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Sound Speed Profile 




Sound Speed (meters/second) 

Figure 3.2; Sound speed profile at receiving array. 



Mode One 



Mode Two 





Eigenfunction Amplitude 

Figure 3.3: Normal modes one and two. 
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Mode Three 



Mode Four 




Eigenfunction Amplitude 

Figure 3.4: Normal modes three and four. 
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Eigenfunction Amplitude 

Figure 3.5: Normal modes five and six. 
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C. DATA ACQUISITION AND PREPROCESSING 

Once deployed, the receiving array recorded numerous data 
sets, each lasting in excess of 60 minutes. Each recording 
comprises 32 channels of acoustic data. The hydrophone output 
was ported to a patch panel and preamplifier where a fixed 
gain of 26 dB was applied. Following this, the signal passed 
through a variable gain amplifier which automatically set the 
signal gain to optimize the 10 volt dynamic range of the 
analog to digital converter. This process had the duel effect 
of enhancing the precision of the data and preventing 
saturation of the A/D converters. Subsequent normalization of 
the data set to hydrophone output levels was accomplished by 
an implementation of a program written by Keith Von der Heydt 
of Woods Hole Oceanographic Institution. 

Analog to digital conversion utilized a sampling frequency 
of 228 Hz. The signal was passed through an analog bandpass 
filter prior to conversion. This prevented aliasing of 
frequencies above Nyquist into the data set. Low frequency 
noise components were also attenuated. The cutoff frequencies 
for the bandpass filter were 10 and 80 Hz. Data storage 
utilized both optical disk and magnetic tape. Data 
acquisition and preprocessing of the data sets used in this 
work was conducted by Miller, Chiu, Frogner and others. 
Figure 3.6 illustrates the equipment alignment. [Ref. 8] 
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Figure 3.6: Data acquisition system 
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V. RESULTS AND CONCLUSIONS 



The desired result of this research is the development of 
a software package capable of bearaforming a tomography signal. 
The beamformer should discriminate among the various modal 
components of the signal, thus enabling one to treat the modes 
individually. Due to the inherent nature of vertical arrays, 
the software needs to incorporate information pertaining to 
the time varying array tilt and individual hydrophone depths. 
To this end, the beamformer utilizes a duty cycle. Each such 
cycle commences with an assessment of array geometry. This 
information is used to determine the appropriate steering 
delays and hydrophone weighting coefficients. Once 
established, the beamformer processes 60 seconds of acoustic 
data before repeating the procedure. In this manner, the 
software package implements a quasi-stable virtual array in 
which the virtual aperture remains fixed in space. 

The results presented examine two performance aspects. 
The first pertains to array geometry description. Calculation 
of the steering delays and amplitude weights is detailed. 
Additionally, the acoustic output of the array is compared to 
that of a single hydrophone on the sound channel axis, thus 
allowing for a quantitative evaluation of the ultimate 
performance of the beamformer. 



26 



A. HYDROPHONE DELAY AND WEIGHTING 



Information regarding array geometry was provided by two 
instrument packages. The upper sensor was located 4.0 meters 
above hydrophone number one, the lower sensor 5.0 meters below 
hydrophone number 20. Figure 5.1 illustrates the instrument 
configuration and coordinate system used in determining the 
array’s spatial orientation. 



NORTH 




Upper Instrument : 

NS current speed 
EW current speed 
Temperature 
Pressure 
Tilt 

Lower Instrument: 

Temperature 
Pressure 
Conductivity 
NS tilt 
EWtilt 



Figure 5.1: Array Instrumentation 
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1. The Upper Sensor 

Information provided by the upper instrument includes 
the component of current velocity in both the north-south and 
east-west directions, ambient pressure, temperature and tilt. 
The direction in which the array tilts is determined by the 
current direction. The primary assumptions being that the 
current acting on the array is independent of depth and that 
the direction in which the array tilts is completely 
determined by the current direction. 

In making these assumptions, one considers the primary 
source of current measured by the instrument to be a result of 
ship positioning adjustments. As mentioned previously, the 
R/V Point Sur was required to periodically move the ship in 
order to remain clear of the array. Such movements were 
intended to straighten the tether between the array and the 
towing rig, thus preventing the tether from becoming fouled in 
the ship's propellers. These adjustments resulted in a 
tensioning of the tether. It is this tensioning that is 
considered to have resulted in array movement and consequent 
water flow across the instrument package. Based on this 
assessment, the assumption that tilt direction was determined 
by the flow measured at the upper instrument is consistent 
with the conditions of the experiment. Higher order models 
could be derived for array dynamics, however this research has 
utilized the first order model of a linear array. 
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The next variable pertinent to array geometry is the 
pressure measured by the upper sensor. From this it is 
straight forward to calculate the sensor's depth using the 
following relationship, 

^ _ P - 1 a tmosphei e 
P9 

where P is the pressure (in pascals) measured by the 
instrument, p is the density of seawater and g is the force of 
gravity . 

Figures 5.2 and 5.3 present the data received from the 
upper instrument package during two of the acoustic recording 
periods. The recording beginning 00:05 (GMT) January 27 
coincides with the reception of a m-sequence of length 2047. 
A continuous wave signal was received during the recording 
beginning 15:05 (GMT) January 27. They are representative of 
the most and least favorable conditions encountered during the 
course of the experiment. 
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Array Geometry 
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Figure 5.2: Array geometry for 00:05 (GMT) January 27. 
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Array Geometry 




Figure 5.3: Array geometry 15:05 (GMT) January 27. 
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2 . The Lower Sensor 



The lower sensor was intended to provide a redundant 
source of array positioning data. Specifically/ it measured 
tilt in the north-south and east-west directions, temperature, 
pressure and conductivity. As discussed previously, the lower 
instrument appears to have suffered a casualty which rendered 
the tilt data suspect. 

Figure 5.4 indicates that an event occurred which 
caused the sensor to provide data which is not consistent with 
the assumed array geometry. One scenario explaining the event 
is that the sensor failed shortly after entering the water. 
Another possible explanation is that the array became 
entangled in itself upon descent. 

All array tilt data was stored for subsequent analysis 
ashore. Therefore, the research team had no knowledge of the 



apparent failure 


during 


the 


conduct of 


the 


experiment . 


As 


there was no 


reason 


to 


suspect 


the 


instrument , 


no 


investigation was made 


to 


determine 


the 


source of 


the 



inconsistencies . 
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Lower Sensor Output 




Time (minutes) 



Figure 5.4: Lower sensor output upon deployment. 
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Depth (meters) 



3. Steering Delays 

Time delays must be applied to each of the hydrophones 
in order to correct for array tilt. However, only the 
component of horizontal displacement collinear to the signal 
propagation direction is used to calculate the time delay. 
This requirement mandates apriori knowledge of the propagation 
direction. Prior work has determined that the expected 
direction from which the signals would arrive is 217.0® true 
[Ref. 9]. 

In addition to array geometry and propagation 
direction, one requires a nominal propagation speed to convert 
the horizontal displacements to time delays. The appropriate 
speed to use is the mode group speed, as previously discussed. 
Group speeds for this study were obtained from eigenvalues 
calculated by the normal mode model courtesy of Chiu and Ehret 
[Ref. 5]. The mode group speed of the carrier frequency was 
used to calculate time delays. Figure 5.5 presents the mode 
group speeds as a function of frequency for the lowest five 
modes . 
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Figure 5.5: Mode group speeds. 
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4. Hydrophone Amplitude Weighting 

Application of hydrophone weights utilizes the array 
geometry to determine the depth of each array element. The 
value of the depth dependent function (Z^) is then assigned to 
the appropriate hydrophone output. This manner of weighting 
enhances the array sensitivity to the m^^ normal mode. 

An additional consideration encountered in hydrophone 
weighting is fault tolerance . Vertical acoustic arrays 
operate in a hostile environment, thus individual hydrophones 
may fail. The amplitude weights applied to the array must 
cope with the failure of individual elements. To this end, 
the beamformer applies a weight of zero to elements which have 
been evaluated as unreliable. 

Log books kept by the research team aboard the R/V 
Point Sur indicating suspect hydrophones were verified by post 
processing statistical analysis of the individual elements. 
Several hydrophones were found to have failed during the 
course of the experiment. Specifically elements 18 through 28 
were evaluated as unreliable. As a result of this 
determination, only hydrophones one through 17 were used to 
beamform the target signals. 

B. ACOUSTIC PERFORMANCE 

The preliminary data analysis seeks to quantify the 
beamformer ' s performance with respect to a single hydrophone 
located on the sound channel axis. The power spectrum 
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calculated from the output of hydrophone five was compared to 
the spectra calculated from the beamformed output of channels 
one through 17. The array gain is defined as 



where SNR^ is the signal to noise ratio of the beamformed 
output and SNR^ is the signal to noise ratio of the n**’ 
hydrophone [Ref. 10]. The SNR for each case was calculated by 
integration of the power spectra over a 6.0 Hz bandwidth 
centered on the continuous wave frequency. 

The array gain was calculated for numerous power spectra, 
each representing 60 seconds of acoustic data. Additionally, 
the gain was calculated for the lowest five normal modes 
present. Values for array gain range from a low of 2.5 dB to 
a high of 8.5 dB. The nominal array gain is 6 dB. 

Reasons for the range in values include temporal 
variability in the modal components of the signal, errors 
introduced in the steering delays and errors in amplitude 
weighting. Time delay and amplitude weight errors are the 
direct result of assumptions made regarding array geometry. 
The primary source of geometric error is considered to involve 
the assumption that tilt direction is completely determined by 
current direction. Array tilt bearing may not be aligned well 
with current direction during periods when the direction of 
flow across the sensor is varying. Therefore, during periods 




( 5 . 2 ) 
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of variable current flow, marginal performance can be expected 
of the beamformer. 

Figures 5.6 through 5.11 present representative power 
spectra of single channel data and beamformed data for modes 
one through five over the entire range of frequencies 
analyzed. 




O 20 40 60 80 100 120 

Frequency (Hz) 

Figure 5.6; Hydrophone five power spectrum. 
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Figure 5.8: Mode two power spectrum. 
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re 5.9: Mode three power spectrum. 
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Figure 5.10: Mode four power spectrum. 
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Mode Five Power Spectrum 
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Figure 5.11: Node five power spectrum. 



As shown in the previous figures, the noise field at the 
beamformer is not incoherent. This fact, coupled with 
imperfect geometric array description, detracts from the 
calculated array gain. 

One feature worth noting is the decreasing strength of the 
undesirable tonal components with increasing mode number, 
specifically those in the 40 to 50 Hz range. Presumably, 
these components were radiated by the R/V Point Sur. One 
hypothesis is that the strength of these signals in mode one 
is the result of grating lobes in the array's beam pattern for 
this mode. The amplitude weights applied to the hydrophones 
are similar to those for a plane wave beamformer, in that 
there is no phase shift applied along the array in the form of 
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coefficients of differing sign. Therefore, spatial aliasing 
is possible for frequencies with wavelengths less than twice 
the hydrophone spacing. For the Heard Island array, this 
includes frequencies ranging from 17 Hz to the bandpass cutoff 
frequency of 80 Hz. 

As the mode number sought increases, these frequency 
components diminish in strength. This indicates that the 
rejection of coherent noise of local origin improves as the 
number of 180° phase changes applied to the hydrophone output 
increases . 

Also shown in these figures is the presence of distant 
sources of coherent noise. These other tonal components are 
seen to vary both as a function of mode number and time. This 
observed behavior is consistent with modal propagation in a 
range dependent environment. Merchant traffic in shallow 
water and offshore drilling activity are two possible 
candidates for sources of coherent noise which becomes coupled 
into the sound channel. 

Figures 5.12 through 5.17 present the same power spectra 
in the immediate vicinity of the carrier frequency (57 Hz). 
Again, the beamformed output indicates successful attenuation 
of undesirable components of a coherent noise field, 
presumably of local origin. 
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Frequency (Hz) 

Figure 5.12: Hydrophone five power spectrum. 
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Figure 5.13: Mode one power spectrum. 
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Frequency (Hz) 

Figure 5.14: Mode two power spectrum. 




Frequency (Hz) 

Figure 5.15: Mode three power spectrum. 
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Mode Four Power Spectrum 




Frequency (Hz) 

Figure 5.16: Mode four power spectrum. 




Frequency (Hz) 

Figure 5.17: Mode five power spectrum. 
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C. CONCLUSIONS 



Modal beamforming is a valuable technique for use in long 
range tomography experiments. The signals presented in this 
work traversed the longest transmission path for which 
reception was attempted during the course of the Heard Island 
Experiment. Although no travel time estimates have been made, 
the successful detection of the signals indicates that large 
scale tomography experiments are feasible. 

The primary concern of this work was the development of 
a system capable of detecting the acoustic signals emanating 
from the vicinity of Heard Island. This has been 
accomplished. Additionally, the data describing the complex 
dynamics of the near vertical acoustic array has implications 
lor the design and construction of vertical arrays for use in 
future tomography experiments. The most troublesome aspect of 
this work has been the accurate determination of array 
geometry, specifically tilt direction. Describing tilt 
direction in terms of measured current at one point on the 
array is considered inadequate. 

D. RECOMMENDATIONS 

Improvements to this implementation include physical 
aspects of the array itself and performance of the numerical 
algorithm. Large vertical arrays are particularly sensitive 
to errors in estimated geometry. In this experiment, an error 
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of 0.25° in tilt (assuming perfect knowledge of tilt 
direction) creates a phase error of 45 degrees between the 
most distant elements. This error is more destructive to the 
higher modes owing to their greater spatial extent. The 
problem is aggravated by questionable tilt direction data. 

Construction of future vertical arrays should incorporate 
multiple instrument packages which directly measure tilt and 
direction. As this was the case with the lower instrument on 
the Heard Island array, it is most unfortunate that it failed. 
Additionally, an investigation into the utility of higher 
order models describing array shape would be illuminating. 
Incorporating a quadratic model for array shape would not 
impose a significant computational load on the beamformer and 
may improve the estimated position of individual hydrophones. 

The current implementation of this beamformer utilizes 
third order polynomial interpolation during the application of 
time delays to individual hydrophone outputs. As such it 
represents the greatest computational load in the program. 
The interpolator is based on Neville's algorithm [Ref. 11]. 
No attempt was made to minimize execution time. As a result, 
the software requires approximately four hours to process one 
hour of acoustic data. A more efficient interpolator may 
reduce execution time sufficiently to permit this technique to 
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be used for real time modal beamforming. To this end, the 
software is assembled from modular components, thus allowing 
for improvements on a function by function basis. 
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APPENDIX A 



The following programs were developed during the course of 
this research. They include the modal beamformer and 
associated utility programs. These programs may be obtained 
by contacting James H. Miller at his address in the initial 
distribution statement. 

A. THE MODAL BEAMFORMER 

1. Operational Considerations 

The modal beamformer was designed with portability a 
prime consideration. As a result, there are numerous program 
parameters which must be set by the end user. Among these are 
the characteristics of the physical array and data acquisition 
system. All such user selectable parameters are contained 
solely within the program definitions. The parameters appear 
below in the same order in which they appear in the program. 

• UNIX VERSION: selectively compiles either unix or ansi 

compatible code blocks. 

• ASCII: selects output in either ascii or binary format. 
Ascii is useful for data export, binary saves storage 
space . 

• SIGNAL: permits user to enable or disable the beamformer 's 
time series output. 

• SPECTRUM: permits the user to enable or disable the 

beamformer 's power spectrum output. 
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• LOWER_SENSOR: directs the program to load or ignore lower 
tilt sensor data. This option should remain OFF as there 
is no reliable lower sensor data. 

• VALIDATE: enables program validation. If selected the 

program will create output files which list the array 
geometry, steering delays, hydrophone weights and 
calculated mode group speed. This feature is useful for 
program debugging. 

• ERROR_ESTIMATE : if selected, this feature causes the 

polynomial interpolator to store the upper bound on the 
mean squared error in the steering delays. This feature 
doubles the output of the beamformer. 

• ON: logical switch for program control. 

• OFF: logical switch for program control. 

• INTERPOLATE: selects interpolator. This feature permits 
the installation of an improved interpolator without the 
requirement to search for and replace each call to the 
function . 

• ORDER: indicates the degree of the polynomial used in the 
polynomial interpolator. 

• STEP: indicates the number of points on either side of a 
sequence to be taken when estimating derivatives. 

• TINY: prevents division by zero in floating point 

operations . 

• PI: used for trigonometric recursions. 

• RADIAN: used for degree to radian conversions. 

• OFFSET: indicates the distance (in meters) between the 

number one hydrophone and the upper tilt sensor. 

• DELTA_R: indicates array element spacing (in meters). 

• CTD_OFFSET: indicates the difference in depth between the 
sound speed profile's first data point and the depth 
increment used in the profile (in meters). 

• SSP_LENGTH: indicates the maximum number of data points 
that an eigenfunction may contain. The only restriction 
applicable to this parameter is available memory. 
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• EIGVAL_LENGTH : indicates the maximurc number of 

eigenvalues that the program will accept. The only 
restriction applicable to this parameter is available 
memory. 

• LOOK_DIRECTION : the compass direction (in degrees true) 

from which the signal arrives. 

• TILT_BUFFER: the buffer length allocated to tilt data. 

The only restriction applicable to this parameter is 
available memory. 

• BUFFER_TIME: the length of time (in seconds) represented 
by the acoustic data input buffer. This period represents 
the output buffer length plus one second on either end to 
permit steering delays and/or advances. The array is 
currently capable of end firing a 1500 meter acoustic 
array. This parameter must be an integer. Additionally, 
60 must be an integer multiple of ( BUFFER_TIME-2 ) . The 
array has a duty cycle of one minute between design 
changes. There must be an integer number of output 
buffers per duty cycle. 

• F_SAMPLE : sampling rate of the data acguisition system. 

This value must be an integer. 

• CHANNELS: the number of hydrophones on the array. 

• F_CARRIER: the carrier frequency of the target signal. 
The group speed used in calculating steering delays is 
based on this frequency. This must be entered in floating 
point . 

• FFT_LENGTH : the length of the FFT used if frequency 

domain output is selected. This value must be a power of 
two which is less than or equal to the number of samples 
in a single channel of input data. 

• SWAP: a macro used in calculating a FFT. 



The following program functions are adaptations of 
those found in Reference 11: 



• polint: 

• realft: 

• fourl: 



performs the polynomial 
calculates Fast Fourier 
calculates Fast Fourier 



interpolation . 
Transforms . 
Transforms 
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• vector: allocates memory for one dimensional arrays. 

• matrix: allocates memory for two dimensional arrays. 

• free_matrix: deallocates memory from two dimensional 

arrays . 

• ExitOnError: provides abnormal program termination. 



2 . Beeunf ormer Source Code 

/irkirk 

* PROGRAM: BEAMFORMER vsn 1.0 

* WRITTEN BY: Steven Crocker 

* LAST UPDATE: October 9, 1991 

* 

★ This program takes input from various data files and the user. It 

* outputs a data file. The inputs are a number of channels of digital 

* acoustic data, and information regarding the physical characteristics 

★ and geometry of the receiving array. Additionally, environmental data 



* in the form of normal mode eigenfunctions and eigenvalues at the 

★ receiving array are required to operate this beamformer. The output 

* is a single channel of acoustic data. 

★ 

Irkirk/ 

^define UNIX VERSION /* either ANSI or UNIX ★/ 

^^define ASCII ON /* select output mode ON or OFF ★/ 

^define BINARY OFF /★ select output mode ON or OFF */ 

^define SIGNAL ON /* either ON or OFF ★/ 

^define SPECTRUM ON /★ either ON or OFF ★/ 

#define L0WER_SENS0R OFF /★ either ON or OFF ★/ 

^define VALIDATE OFF /* either ON or OFF */ 

#define ERROR^ESTIMATE OFF /★ either ON or OFF ★/ 

^define ON 1 /* logical "switch" */ 

^define OFF 0 /* logical "switch" */ 

^define INTERPOLATE polint /* polint */ 

^define ORDER 3 /* order of interpolator (odd)*/ 

//define STEP 1 /* number of steps for derivatives */ 

//define TINY 1.0e-25 /* prevents division by zero . */ 

//define PI 3.14159265359 /* for freq to omega conversions */ 

#define RADIAN 57.2957795131 /* for degree to radian conversions */ 

//define OFFSET 4.0 /* dist btwn upper sensor and phone //I */ 

//define DELTA_R 45.0 /* array element spacing */ 

//define CTD^OFFSET 0.5 /* diff btwn ctd depth inc & 1st depth */ 

#define SSP_LENGTH 2500 /* max number of pts in eigunfunction */ 

//define EIGVAL_LENGTH 230 /* max number of eigenvalues */ 

//define LOOK_DIRECTION 217.0 /* direction from which signal arrives */ 

//define TILT_BUFFER 120 /* max length of tilt data vectors */ 

#define BUFFER_TIME 12 /* input buffer length in seconds (int)*/ 

//define F_SAMPLE 228 /* sampling frequency (int)*/ 

//define CHANNELS 32 /* number of channels processed */ 

//define F_CARRIER 57.0 /* carrier frequency */ 

#define FFT_LENGTH 2048 /* radix 2 <= (BUFFER_TIME-2)*F_SAMPLE */ 

//define SWAP(a,b) tempr=(a); (a) = (b); (b)=tempr 



//include<stdio. h> 
//include<mal loc. h> 
//include<math.h> 

#if defined ( ANSI ) 
//include<f loat . h> 
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#incLude<stdlib.h> 
int getlnput(void); 
int putOutput (void); 

int processTi It (f loat float float 

int processModesC float float *weight, float *ptrC); 

int dydxCfloat *x, float *y, float *ddx, int points); 

int polintCfloat *xa, float *ya, int n, float x, float *y, float *dy); 

int realftCfloat *data, int n, int isign); 

int fourKfloat *data, int nn, int isign); 

float *vector(int length); 

float **inatrix(int row, int col); 

int free_matrix (float int row); 

void ExitOnError(char error_txtC]); 

#elif defined ( UNIX ) 
int vectorO, 

■latrixO; 

get Input (), * 

putOutput (), 
processTi It (), 
processModesO, 
dydxO, 
polintO, 
realftO, 
fourlO, 
f ree^roatrixO, 

ExitOnError(); 

#endif 

/* Global Variable Declarations */ 
float **inSOUND, *outSOUND, *MeanSqError, Max=0.0, Min=1.0e25; 
int lastMinute, Minute=1, f irstBuffer=1; 

mainO 

C 

FILE *fpVl; 

int i, j, k, n, *shift; 

float *^x, **y, *indx, *satnples, arg, ans, err, *delay, 

★delta, *weight, *pwrSpectrum, Cgroup; 



/* Memory Allocation and Tilt Data Processing ★/ 
x=(f loat**)matrix(CHANNELS, TILT^BUFFER); 
y=(f loat**)matrix(CHANNELS, TILT^UFFER); 
z=(f loat**)matrix(CHANNELS, TILT__BUFFER); 
processTi It (x, y, z); 

if (VALIDATE==ON) 

printf ("\nDumping array geometry to array .dat\n"); 
f pV1=f open ("a rray.dat ", "wt"); 

if (fpV1==NULL) ExitOnErrorC'Error opening validation file"); 

fprintf (fpV1 ,"Channel\t\t X\t\t Y\t\t Z\n"); 

for (i=1; i<=lastMinute;i++) 

fprintf(fpV1, "MINUTE: %i\n", i); 
for (j=1; j<=CHANNELS;j++) 

fprintf(fpV1,"%i\t %f\t %f\t %f\n", j,xCj]Ci],yCj]Ci3,zCj]Ci]); 
> /* for */ 
f close(fpVl); 

> /* if */ 

/★ Memory Allocation*/ 
weight=(f loat*) vector (CHANNELS); 



53 



indx=(f Loat*)vector(OROER+1 ); 
delay=(f Loat*)vector(CHANNELS); 
delta=(f loat*)vecto^(CHANNELS); 

outSOUND=(f Loat*)vecto^((BUFFER_TIME-1)*F_SAMPLE); 

^nS0UND=(f Loat**)matrix(CHANNELS, BUFFER_TIME*F_SAMPLE); 
if ( (shif t=(int*)maUoc( (CHANNELS+1)*si2eof (int)) )==NULL) 

Exi tOnErrorC'Memory allocation failure for shiftC]."); 
for (i=1; i<=ORDER+1; i++) indxCi]=(f loat)i; 
if (ERROR EST1MATE==0N) 

MeanSqError=(f loat*)vector((BUFFER_TIME-1)*F_SAMPLE); 
if (SPECTRUM==ON) pwrSpectrum=(f loat*)vector(FFT_LENGTH/2); 
whi le(Minute<=lastMinute) 



/* Mode Data Processing */ 
processModesCz, weight, &Cgroup); 

/* Calculate delays, shifts, etc */ 
for (i=1; i<=CHANNELS; i++) 

C 

delayCi ]=x[i][Minute]/Cgroup; /* Time delay */ 

shiftCi]=(int)(delay[i]*(f loat)F_SAMPLE); /* # of samples */ 

/* fraction of 1 sample */ 
deltaCi]=delay[i]*(f loat) F_SAMPLE-( float )shift[i]; 

> /* for */ 

if (VALIDATE==ON) 

printf ("Dumping element delays to delay. dat\n"); 

if (f i rstBuf fer) 

C 

if ( (fpVl=fopen("delay.dat", "wt")) == NULL) 

ExitOnErrorC'Error opening validation file."); 

> /* if */ 
else 
{ 

if ((fpVl=fopen("delay.dat", "at")) == NULL) 

ExitOnErrorC'Error opening validation file."); 

> /* else */ 

f printf (fpVl , "MINUTE: %i\n". Minute); 
fprintf (fpVl , 

"ChannelXt delay\t\t int shiftXt fraction of 1 shiftXn"); 

for (i=1;i<=CHANNELS;i++) 

fprintf (fpVl, "%i\t %e\t %i\t% e\n", 
i ,delayCi],shiftCi3,deltaCi]); 

> /* for */ 

f c lose(f pvl ) ; 

printf ("Dumping phone weights and grp speed to modal .dat\n"); 

if (f i rstBuffer) 

{ 

if ( (fpVl=fopen("modal.dat", "wt")) == NULL) 

ExitOnErrorC'Error opening validation file.Xn"); 
fprintf (fpVl , "Group speed for F_CARR1ER is: %g\n\n\n", Cgroup); 

> /* if */ 
else 

if ((fpV1=fopen("modal .dat","at")) == NULL) 

ExitOnErrorC'Error opening validation file.Xn"); 
fprintf (fpVI , "XnHydrophone weights for minute %i\n", Minute); 

> /* else if */ 
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fprintf (fpVl, "Channel \tWeight\n"); 

for (i=1;i<=CHANNELS;i++) 

fprintf (fpVI, "%i\t%e\n", i,we^ghtCi]); 

f closeCf pVI); 

> /★ if ★/ 

if (SPECTRUM==ON) 

for (k=1;k<=FFT_^LENGTH;k++) pwrSpectrutnCk]=0.0; 

for (n=1; n<=60/(BUFFER TIME-2); n++) 

C 

getInputO; 

if (f i rstBuf f er) /* Produce first output buffer ★/ 

T 

for (i=1; i<=F_SAMPLE; i++) 

i 

outSOUNDCi]=0.0; 

if (ERROR_^ESTIMATE==ON) MeanSqErrorCi ]=0.0; 

> /* for */ 

for (i = F_^SAMPLE+1; i<=(BUFFER_TIME-1)*F_SAMPLE; i++) 

i 

outSOUNDCi]=0.0; 

if (ERROR_ESTIMATE==ON) MeanSqErrorCi ]=0. 0; 

for (j=1; j<=CHANNELS; j++) 

i 

if (delayCj]>=0.0) 

C 

arg=(f loat) (ORDER+1)/2.0+(1-deltaC j]); 
samples = &inSOUNDCj][i-shiftC j]-(0RDER+1)/2]; 

> /★ if ★/ 

else if (delayCj]<0.0) 

< 

arg=(f loat) (0RDER+1)/2.0+delta[j]; 

samples = &inS0UNDCj][i-shiftCj]-(0RDER+1)/2+1]; 

> /★ else if ★/ 

INTERPOLATECindx, samples, ORDER+1, arg, &ans, &err); 
out S0UNDCi]=out SOUND Ci]+ans; 
if (ERROR_^ESTIMATE==ON) 

MeanSqErrorCi]=MeanSqErrorCi]+err*err; 

> /★ for */ 

if (fabs(outSOUNDCi])>Max) Max=fabs(outSOUNDCi]); 
if (fabs(outSOUNDCi])<Min) Min=fabs(outSOUND[i]); 
if (ERR0R_^ESTIMATE==0N) 

MeanSqErrorCi]=MeanSqError[i]/(f loat ) CHANNELS; 

> /★ for ★/ 

> /* If ★/ 

else /* Produce subsequent output buffers */ 

C 

for (i = F_^SAMPLE+1; i<=(BUFFER_TIME-1)^F__SAMPLE; i++) 

i 

outS0UNDCi-F_^SAMPLE]=0.0; 

if (ERROR_ESTIMATE==ON) MeanSqErrorCi-F_SAMPLE]=0.0; 

for (j=1; j<=CHANNELS; j++) 

if (delayCj]>=0.0) 

C 

arg=(f loat) (ORDER+1) /2.0+(1-deltaCj]); 
samples = &inSOUNDCj][i-shif tCj]-(0RDER+1)/2]; 

> /★ if ★/ 

else if (delayCj]<0.0) 

C 

arg=(f loat)(ORDER+1)/2.0+deltaCj]; 

samples = &inS0UNDCj][i-shiftCj]-(0RDER+1)/2+1]; 
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y /* else if */ 

INTERPOLATECindx, samples, ORDER+1, arg, &ans, &err); 
outSOUND[i-F_SAMPLE]=outSOUNDCi-F_SAMPLE]+ans; 
if (ERROR_ESTIMATE==ON) 

MeanSqErrorCi-F_SAMPLE]=MeanSqErrorCi-F_SAMPLE]+err*err; 
> /* for */ 

if (fabs(outSOUND[i-F_SAMPLE])>Max) 
Max=fabs(outSOUNDCi-F_SAMPLE])/ 

if (fabs(outSOUNDCi-F_SAMPLE])<Min) 
Min=fabs(outSOUND[i-F_SAMPLE]); 
if (ERROR_ESTIMATE==ON) 

MeanSqE r ror [ i -F_SAMPLE]=MeanSqEr ror [ i -F_SAMPLE3 / 

(float) CHANNELS; 

> /* for */ 

> /* else */ 

if (SIGNAL==ON) putOutputO; 

if (SPECTRUM==ON) 

C 

windowCoutSOUND, FFT_LENGTH); 
realftCoutSOUND, FFT_LENGTH/2, 1); 
for (k=0;k<FFT LENGTH;k += 2) 

pwrSpectrumCk/2+13=pwrSpectrumCk/2+13+ 

outSOUND[k3*outSOUNDCk3+outSOUNDCk+13*outSOUND[k+13; 

> /* for */ 

y /* if */ 

f i rstBuffer=0; /* set firstBuffer to false */ 



> /* for */ 

if (SPECTRUM==ON) 

for (k=1; k<=FFT^LENGTH; k++) 

pwrSpectrum[k3=pwrSpectrumCk3*(f loat ) (BUFFER_TIME-2)/60.0 

dumpSpectrum(pwrSpectrum); 

> /★ If */ 

printf("\t%i minutes of input data processed. \n". Minute); 
Minute++; /* increment minute counter */ 

> /* whi Le */ 

printf ("EXECUTION COMPLETE: End of tilt data encounteredXn"); 
printf ("Maximum magnitude encountered was: %e\n",Max); 
printf ("Minumum magnitude encountered was: %e\n",Min); 
exit(O); 

y /'kic’k’kticicirkirkici^irki^irki^'kiric-k'ki^'k-k'k'k'k'k'k'k'k'kirk'k-kirkic'k'k'k’k-k-k'k t 'k'k * * k 'ki r kirki r ki r kt'k f 
f 'kiricic'k'kicicicic'k iciiicic'k'k'kirkirk'k'k'k'k END main irkiciri^i^i^ticifi^iririri L 'k- k 'k'k'k'k-k'k'k'k'k'k'k'k'k / 

/irkirk^ 

★ FUNCTION: getInputO 

* 

* This function handles all acoustic input. It also provides one of 

* two normal process terminations available in the program. (The other 

* is located in mainO.) 

* 

* Arguments: none 



* Return value: 

★ 

* Functions called: 

* 

* Definitions called: 



0 



vector () 



ANSI 

F_SAMPLE 
BUFFER TIME 



UNIX 

CHANNELS 



ExitOnErrorO 



* 

★ 

* 



* Global variables called; 



inSOUNDC3[3 

firstBuffer 



Min 

Max 



★ 

* 
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* Significant memory allocation: diskBufferC] 

★ 

tirkirk/ 

ttM defined ( ANSI ) 
int getlnput(void) 

#elif defined ( UNIX ) 
get Input () 

#endif 

int i, j, buffer, items; 
float *-diskBuffer; 
char fileNameCSOD; 
static FILE *^fpStatic; 

if (f i rstBuffer) 

-C 

printf ("Enter file name for input acoustic data: "); 

if((scanf("*^s", fi leName) )==EOF) 

Exi tOnErrorC'Fatal error in scanfO"); 

printf ("\n\n"); 

if ( (fpStatic=fopen(fi leName, "rb")) == NULL) 

ExitOnErrorC'Error opening INPUT ACOUSTIC data file."); 
> /★ if */ 



/* Memory Allocation */ 

if (firstBuffer) buf fer=BUFFER TIME*F_SAMPLE*CHANNELS; 
else buffer=(BUFFER_TIME-2)*rF“SAMPLE*CHANNELS; 
diskBuf fer=(f loa t*) vector (buffer); 



items=f read( (char*) (diskBuffer+1), sizeof (float), buffer, fpStatic); 
if (items==buffer) /^continue*/; 
else if (f error(fpStati c) != 0) 

ExitOnErrorC'Error encountered while reading input acoustic data"); 
else if (feof (fpStatic) != 0) 



pr i ntf ( " \ n \ t ************************************************** \ n " ) ; 
printf ("\t End of File reached: EXECUTION COMPLETEXn" ); 
printf ("\t\t%i minutes of data processedXn", Minute-1 ); 
printf ("\t\t%i bytes of data discardedXn", items*si zeof (f loat) ); 
printf ("\t\tMaximum magnitude encountered was: %e\n",Max); 
printf ("\t\tMinimum magnitude encountered was: %e\n",Min); 
printf ("\t End of File reached: EXECUTION COMPLETEXn"); 

fclose(fpStatic); 

exit(O); 

> /* else if */ 

else Exi tOnErrorC'Unknown error handling acoustic input file."); 



for (i=1; i<=BUFFER_TIME*F SAMPLE; i++) 

i 

for (j=1; j<=CHANNELS; j++) 

if (firstBuffer) 

{ 

inSOUND[j]Ci] = diskBufferCCHANNELS*(i-1)+j3; 

> /* if */ 
else 

£ 

if(i<=2*F_SAMPLE) 

inS0UNDCj]Ci]=inS0UNDCj][i+(BUFFER_TIME-2)*F_SAMPLE]; 

else 

inS0UNDCj]Ci]=diskBufferCCHANNELS*(i-2*F_SAMPLE-1)+j]; 

> /* else */ 

> /* for */ 

> /* for */ 
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> 



/* Deallocate Memory */ 
f ree( (cha^*)d^skBuffe^); 
^etu^n( 0 ); 



/★★★★★ end getinput ’ktirtrk/ 



* FUNCTION: putOutputO 



This function handles all acoustic output. Additionally, it outputs 
the estimated mean squared error from the interpolators (if enabled). 



* 

* 

* 

* 

★ 

★ 

* 

* 

* 

* 

★ 

* 

* 

* 

* 

* 

★★★★* / 

#if defined ( ANSI ) 
int putOutput(void) 
#elif defined ( UNIX ) 



Arguments: 

Return value: 

Functions called: 
Definitions called: 

Global variables called: 



none 

0 

ExitOnErrorO 

ANSI 

F_SAMPLE 

outSOUNDED 
f irstBuf fer 



Significant memory allocation: none 



UNIX 

BUFFER_TIME 

MeanSqErrorCD 



putOutputO 

#endif 

int i, cut; 

char f i leName[123, mode[23; 
static FILE *fpOutSound, *fpMSE; 



if (f i rstBuf fer) 

if (ASCII==ON) 

mode[03=' W ; 
modeCI 3= ' t ' ; 

> /* if */ 

else if (BINARY==ON) 

mode[03=' W; 
modeC13=' b' ; 

> /* else if */ 



cut=1; 

printf ("Enter file name for output data: "); 

if ( (scanf ("%s", f i leName) )==EOF) 

Exi tOnErrorC'Fatal error in scanfO"); 

printf ("\n\n"); 

if ( (fpOutSound=fopen(fi leName, mode)) == NULL) 
ExitOnErrorCError opening OUTPUT data file."); 

if (ERROR^ESTIMATE==ON) 

{ 

printf ("Opening file error .dat\n\n"); 

if ((fpMSE=fopen("error.dat", mode)) == NULL) 
ExitOnError ("Error opening error.dat"); 

> /* if */ 
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> /★ if */ 
else cut=2; 



if (ERROR_ESTIMATE==ON) 

if (ASCII==ON) 
i 

for (i=1; i<=(BUFFER_TIME-cut)*F_SAMPLE; i++) 
fprintf (fpMSE, "%e\n", MeanSqErrorCi]); 

> /★ if */ 

else if (B1NARY==0N) 
i 

if (fwri te((char*)(MeanSqError+1),si2eof (f loat), (BUFFER TlME-cut )* 
F_SAMPLE,fpMSE)==(unsigned)(BUFFER^TlME-cut)*F_SAMPLE) ; 
else if (ferror(fpMSE) != 0) 

ExitOnErrorC'Error encountered writing error data"); 
else ExitOnError ("Unknown error handling error file."); 

> /★ else if ★/ 

> /* if */ 



if (ASC11==0N) 

< 

for (i=1;i<=(BUFFER_TIME-cut)*F_SAMPLE; i++) 
f printf (f pOutSound, "%e\n", outSOUNDCi]); 

> /★ if ★/ 

else if (B1NARY==0N) 

C 

if (f write ((char*) (out SOUND+1), si zeof (float), (BUFFER_TIME-cut)* 
F_SAMPLE,fpOutSound)==(unsigned)(BUFFER_TIME-cut)*F_SAMPLE) ; 
else if (f error(fpOutSound) 1= 0) 

ExitOnErrorC'Error encountered writing output acoustic data"); 
else Exi tOnErrorC'Unknown error handling acoustic output file."); 

> /* else if */ 
return( 0 ); 

y /iilt'k'k'k'kk'k'k'k'k'k'k'k'k'kirk'k'kic-k'k'k-k j 
firkk-kk ENQ putOutpUt ***★*/ 

/irkirkkirkirkirk*kirki:i:ki^^ / 

/■kirkkk 



* FUNCTION: dumpSpect rumO 

★ 

* This function handles dumps the 

★ 

* Arguments: 

★ 

* Return value: 

★ 

* Functions called: 

★ 

* Definitions called: 

★ 

★ 

* Global variables called: 

* 

* Significant memory allocation: 

★ 

'kkkt'k / 

#if defined ( ANSI ) 
int dumpSpectrum ( float *pwrSpect 
#elif defined ( UNIX ) 
dumpSpectrum ( pwrSpectrum ) 
float *pwrSpectrum; 

#endif 

i 



signal power spectrum (if selected). 
pwrSpectrum 
0 

Exi tOnErrorO 

ANSI UNIX 

F_SAMPLE FFT_LEN6TH 

none 

none 

im ) 



int i; 

static float sequence=1 .0; 
char fileName[12], mode[2]; 
static FILE *fp; 
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if (f i rstBuffer) 

if (ASCII==ON) 

modeCO]=' W ; 
modeC1]='t *; 

> /★ if ★/ 

else if (BINARY==ON) 

modeCO]='w’; 

modeC1]=‘b'; 

> /* else if */ 

printf ("Enter file name for output SPECTRUM data: "); 

i f ( ( scanf ( "%s " , f i leName) )==EOF ) 

ExitOnErrorC'Fatal error in scanfO"); 
printf ("\n\n"); 

if ( (fp=fopen(fi leName, mode)) == NULL) 

ExitOnErrorC'Error opening OUTPUT SPECTRUM data file."); 

printf ("Select output format: \n"); 
printf (" Enter 0 for MATLAB compatible output. \n"); 
printfC Enter 1 for GRAFTOOL compatible output.Xn"); 
if((scanf("%i", &i))==E0F) 

ExitOnErrorC'Fatal error in scanfO"); 
if((i != 0) && (i i= 1)) Exi tOnErrorC'Invalid output selection") 
> /* if */ 



if (ASCII==0N) 

if (fi rstBuffer format==1) 

•C 

fprintf (fp,"0 "); 

for (i=1, i<=FFT LENGTH/2; i++) 

fprintf (fp,"%e ", (f loat) ( i-1 )*F_SAMPLE/FFT_LENGTH); 
fprintf (f p, "\n"); 

> /* if */ 



if (format==1) fprintf (fp, "%e ", sequence); 

for (i=1;i<=FFT_LENGTH/2; i++) 

fprintf (fp,"%e ", pwrSpectrumCi]); 



> 



fprintf (fp,"\n "); 
$equence=sequence+1 .0; 
> /* if */ 

else if (BINARY==0N) 

{ 



if (fwrite( (char*) (pwrSpectrum+1), si zeof (float), FFT_LENGTH/2,fp)= 
(unsigned)FFT_LENGTH/2) ; 
else if (f error(fpOutSound) != 0) 

ExitOnErrorC'Error encountered writing output SPECTRUM data") 
else Exi tOnError ("Unknown error handling acoustic SPECTRUM file."); 

> /* else if */ 
return( 0 ); 



/irk^irk end dumpSpectrum **★**/ 



/irkirkk 

* FUNCTION: processTi It ( ) 

★ 

* This function handles all array tilt data. It calculates the X, Y, Z 

* coordinates of each hydrophone as a function of time. The coordinate 

* system is oriented such that X points toward the signal "origin" and 

* Z points down. 

* 
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* Arguments: 



xC][] 

zC][] 



yC][] 



★ 

* 



* Return value: 

* 

* Functions called: 



0 



Exi tOnErrorO matrixO 
vectorO free matrixO 



★ 

★ 



* Definitions called: 



ANSI 

DELTA_R 

RADIAN 



UNIX 

LOWER_SENSOR 

CHANNELS 



★ 

★ 

★ 

★ 

★ 



LOOK^DIRECTION OFFSET 
TILT BUFFER 



★ Global variables called: 

★ 

* Significant memory allocation: 



lastMinute 



★ 

★ 

★ 

•k 



XX[][] yyC][] 

zzC][] tiltC] 

angleC] udepthC] 

IdepthC] 



kkkkk ! 

#if defined ( ANSI ) 

int processTi lt(f loat *^x, float **y, float **z) 

#elif defined ( UNIX ) 
processTi lt(x,y,z) 
float **x, **y, **z; 

^endif 

int i, j, not EOF; 

float *tilt, ★angle, *udepth, ★Idepth, **xx, **yy, **zz, theta; 
char fileNameC12]; 

FILE *fp1, ★fp2; 

/★ Open Data Files ★/ 

printf ("Enter file name for upper tilt sensor data: "); 

if((scanf("%s", f i leName) )==EOF) 

ExitOnError ("Fatal error in scanfO"); 

printf ("\n\n"); 

if ( (fp1=fopen(fi leName, "rt")) == NULL) 

ExitOnErrorC'Error opening UPPER TILT data file."); 

if (LOWER SENS0R==0N) 

printf ("Enter file name for lower tilt sensor data: "); 

if ((scanf ("%s", f i leName) )==EOF) 

ExitOnErrorC’Fatal error in scanfO"); 

printf ("\n\n"); 

if ((fp2=fopen(fi leName, "rt")) == NULL) 

ExitOnErrorC'Error opening LOWER TILT data file."); 
ldepth=(f loat*)vector(TILT BUFFER) ; 

> /★ if ★/ 

/★ Memory Allocation ★/ 
ti lt=(float*)vector(TILT_BUFFER); 
angle=(float^)vector(TILT_BUFFER); 
udepth=( float*) vector (TILT BUFFER); 
xx=(float**)matrix(CHANNELS, TILT_BUFFER); 
yy=(f loat**)matrix(CHANNELS, TILT^BUFFER); 
zz=(f loat**)matrix(CHANNELS, TILT“buFFER) ; 

i=1; /* read upper tilt data ★/ 
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notE0F=1; 
whi le(notEOF) 

if (fscanf (fp1/‘%g %g %g\n",&ti LtCi!l,&angLe[i3,8iudepth[i]) != EOF) i++; 
else notE0F=0; 

> 



if (LOWER_SENSOR==ON) 

i 

j=1; /* read lower tilt data */ 

notE0F=1; 
whi le(notEOF) 
i 

if (fscanf(fp2, "%g\n",&ldepth[j]) != EOF) j++; 
else notE0F=0; 

> /★ while */ 

if (i<=j) lastMinute=i-1; 

else lastMinute=j-1; 

> 

else lastMinute=i-1; 



/************This is the assumed array geometry: LINEAR***********/ 
for (j=1; j<=CHANNELS; j++) 



for (i=1; i<=lastMinute; i++) 

i 

xx[j][i]=DELTA_R*(f loat) ( j-1 )*cos(ang LeCi D/RADIAN)* 
sinCti ItCiD/RADlAN); 

yyCj ]Ci]=DELTA_R*(f Loat) (j-1 )*sin(angleCi]/RADIAN)* 
sin(tilt[i]/RADlAN); 

zzCj][i]=DELTA_R*(f loat) (j-1 )*cos(ti ltCi]/RADIAN)+ 
OFFSET*cos(tilt[i]/RADlAN)+udepthCi]; 

> /* for */ 

> /* for */ 



theta=(360.0-LOOK_DlRECTION)/RADIAN; /* coordinate rotation */ 
for (j=1; j<=CHANNELS; j++) 

for (i=1; i<=LastMinute; i++) /* points x into signal */ 

xCj]Ci3=xxCj]Ci]*cos(theta)-yyCj][i]*sin(theta); 
yC j]Ci 3=xxCj 3Ci3*sin( theta )+yyCj 3 [i]*cos (theta); 
2[j3[i3=zz[j3[i3; 

> /* for */ 

> /* for */ 



/* Memory Deallocation */ 
if (LOWER_SENSOR==ON) free((char*) (depth); 
f ree((char*)ti It); 



f ree( (char*)angLe); 
f ree( (char*)udepth); 
f ree_matrix(xx, CHANNELS); 
f ree_matrix(yy, CHANNELS); 
f ree_matrix(zz, CHANNELS); 
f close(fp1 ) ; 

if (LOWER_SENSOR==ON) fclose(fp2); 
return( 0 ); 

/*irkirk end processTilt ***★*/ 



/irk-kirk 



* FUNCTION: processModes( ) 

* 

* This function handles the normal mode data. It calculates hydrophone 

* weights and group speed. The user must insure that the depth vector 

* and eigenfunction vector are of equal length. 
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* 

* 

* 

ic 


Arguments: 


zc:c] 

ptrC 


weightC] 


* 


Return value: 


0 




* 


Functions called: 


vector( ) 


Exi tOnError( ) 


* 




INTERPOLATE () 


dydxO 


* 


Definitions called: 


ANSI 


UNIX 


* 




PI 


SSP LENGTH 


* 




EIGVAL_LENGTH 


ORDER 


* 




CHANNELS 


F_CARRIER 


* 


Global variables called: 


Minute 




* 


Significant memory allocation; 


depthCD 


ZnC] 


* 


OnOffC] 


wc: 


* 

* 

irkirirk/ 

#if defined ( ANSI ) 


KrC] 


dwdK 


int 


processModes(f loat **z, float *weight, float *ptrC) 



#elif defined ( UNIX ) 
processModes(z,weight,ptrC) 
float *weight, *ptrC; 

^ end if 

int i, j, ptsEigVal, set, notEOF, deadPhones, weightNotAssigned; 
static int ptsEigFun; 
float *w, *Kr, *dwdK, err; 

static float depth[SSP_LENGTH+1 ], Zn[SSP_LENGTH+1 ], OnOf f CCHANNELS+1 ]; 
char key, f i leNameC12]; 

FILE *fp1, *fp2; 

if (f i rstBuffer==1 ) 

w=(f loat*)vector(EIGVAL LENGTH); 

Kr=(f loat*) vector (EIGVAL_LENGTH), • 
dwdK=(float*)vector(EIGVAL_LENGTH); 

printf ("Enter file name for normal mode data (eigenfunction): "); 

if ((scanf ("%s", f i leName) )==EOF) 

ExitOnErrorC'Fatal error in scanfO"); 

printf("\n\n"); 

if ((fp1=fopen(fi leName, "rt")) == NULL) 

ExitOnErrorC'Error opening NORMAL MODE data file (eigenfunction)"); 

printf ("Enter file name for normal mode data (eigenvalues): "); 

if ((scanf ("%s", f i leName) )==EOF) 

ExitOnErrorC'Fatal error in scanfO"); 

printf("\n\n"); 

if ((fp2=fopen(fi leName, "rt")) == NULL) 

ExitOnErrorC'Error opening NORMAL MODE data file (eigenvalues)."); 

i=1; /* read normal mode data */ 

notE0F=1; 

whi le(notEOF) 

C 

if(fscanf(fp1, "%g %g \n", &depth[i], &ZnCi]) != EOF) i++; 
else notE0F=0; 

> 
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ptsEigFun=i-1; 

for(i=1; i<=ptsEigFun; i++) depth [-i]=depth[ i]+CTD_OFF SET; 
i=1; 

notE0F=1; 
whi le(notEOF) 

if(fscanf(fp2, "%g %g \n", gwCi], &Kr[i]) != EOF) i++; 
else notE0F=0; 

} 

ptsEigVal=i-1; 

for (i=1;i<=CHANNELS;i++) OnOff [i]=1 .0; 

printfCDo you want to turn off any hydrophones? "); 

if ( (scanf ("%s",&key))==E0F) 

ExitOnErrorC'Fatal error in scanfO"); 

if(key=='y’ II key=='Y') 

printf ("\nHow many hydrophones must be secured? "); 

if ((scanf &deadPhones))==E0F) 

ExitOnErrorC'Fatal error in scanfO"); 

for(i=l ; i<=deadPhones; i++) 

printf ("\nEnter hydrophone number to secure; "); 

if((scanf("%i", &j))==E0F) 

ExitOnErrorC'Fatal error in scanfO"); 

if (j>CHANNELS I I j<1) 

Exi tOnErrorC'Bad hydrophone identification"); 

OnOff [j]=0.0; 

> /* for */ 

> /* if */ 

} /* If */ 

for(i=1 ; i<=lastMinute; i++) 

if (depth[ptsEigFun]<z[CHANNEL$][i]) 

printfC'Max eigenfunction depth is; %f \n",depth[ptsEigFun]); 
printfC'at depthC] index of: %i \n",ptsEigFun); 
printfC'Max depth of phone number %i is: %f\n\n", 

CHANNELS, zECHANNELS][i]); 

ExitOnErrorC'Fatal data set error"); 

> /* if */ 

> /★ for */ 

j=i; 

for (i=1; i<=CHANNELS; i++) 
weightNotAssigned=1 ; 

while < j<=ptsEigFun && weightNotAssigned) 

C 

if (z[i][Minute]<0.0 !l depthCj ]<0.0) 

C 

printf("i=%i\n",i); 

printf ("j=%i\n",j ); 

printf("Minute=%i\n",M^nute); 

printf ("z[i][Minute] is: %t\n", zCi][Minute]); 

printf C'depthCj] is: %f \n",depthCj]); 

printf ("Depth less than zero encountered in processModes. "); 
printf ("\n\n"); 
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Exi tOnErrorC'Check input depths for coordinate orientation") 
> /★ if ★/ 



if (depth[j]<zCi][Minute] && depth[j+1]>z[i]CMinute]) 
set=(0RDER+1)/2; 

INTERP0LATE(&depth[j-set3, &ZnCj-set], ORDER+1, z[i][Minute] 
&weight[i], &err); 
weight Not Ass igned=0; 

> /* if ★/ 
j++; 

> /★ while ★/ 

> /* for */ 

for (i=1; i<=CHANNELS; i++) weight[i]=OnOf f Ci]*weight[i]; 

if (Minute==1) 

dydx(Kr,w,dwdK,ptsEigVal); 
for (i=1; i<=ptsEigVal; i++) 

C 

if (w[i3<2.0*PI*F_CARRIER && w[i+1 ]>2 . 0*PI*F_CARRIER) 
set=(0RDER+1)/2; 

INTERPOLATE(&w[i-set],&dwdK[i-set], ORDER+1, 
2.0*PI*F_CARRIER,ptrC,&err); 

> /★ if ★/ 

> /★ for */ 



f ree((char*)w); 
f ree( (char*)Kr); 
f ree((char*)dwdK); 

> /* if ★/ 
returnC 0 ); 

} /irir'k'k'k'k'kic'k'k'k' k 'k k 'k'k'k'k'k'kiric'k'kirk'kic / 



/irkirkif END processModes irk-k^*/ 

/irk kkkirkkkirkkk irkirkk kkkkkkkkkk j 
/krkkirk 



★ FUNCTION: dydxO 

* 

* This function estimates derivatives. 



* Arguments: 

★ 

★ 

★ 

* 
k 
k 
k 
k 
k 
k 
k 
k 
k 

kkkkk/ 

#if defined ( ANSI ) 



Return value: 

Functions called: 

Definitions called: 

Global variables called: 
Significant memory allocation: 



x[] 

ddx[] 



ExitOnErrorO 

ANSI 

STEP 

none 

none 



int dydxCfloat *x, float *y, float *ddx, int points) 

#elif defined ( UNIX ) 

dydx(x,y,ddx, points) 

float *x, *y, *ddx; 

int points; 

#endif 



int n; 



yC3 

points 



UNIX 
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for (n=1;n<=points;n++) 

if ((n>=STEP) gtS (n<=points-STEP)) /*center^/ 

ddx[n3=(yCn+STEP]-y[n-STEP])/(x[n+STEP]-xCn-STEP3); 

else if (n<STEP) /^beginning*/ 

ddx[n]=(y[n+STEP3-yC1])/(xCn+STEP3-x[13); 

else if (n>points-STEP) /*end*/ 

ddxCn3=(y[points3-yCn-STEP3)/(xCpoints3-xCn-STEP3); 

else 

ExitOnError("Index error in dydx"); /★ sanity check */ 
> /* for */ 
returnC 0 ); 

> 

/irkirkic END dydx irkirkk/ 

/kkirkk 

* FUNCTION: poLintO 

* 

* This function performs polynomial interpolation. 



* 

* Arguments: 


xaC3 


yaC3 


* 


n 


X 


* 


y 


dy 


* 

* Return value: 


0 




★ 

* Functions called: 


vectorO 


ExitOnErrorO 


* Definitions called: 


ANSI 


UNIX 


* Global variables called: 


none 




* Significant memory allocation: 

i( 


d[3 


cC3 


★****/ 

#if defined ( ANSI ) 

int polintC float *xa, float *ya. 


int n, float x. 


float *y, float *dy) 



#elif defined ( UNIX ) 
polint(xa,ya,n,x,y,dy) 
float *xa, *ya, x, *y, *dy; 
int n; 

#endif 

C 

int i, m, ns=1; 

float den, dif, dift, ho, hp, w; 
float *c, *d; 

dif=fabs(x-xa[1]); 

c=(f loat*) vector (n), • 
d=(float*)vecto^(n); 

for (i=1; i<=n; i++) 

C 

if ( (dif t=fabs(x-xaCi3) ) < dif) 

ns=i; 

dif=dift; 

> /* if */ 

c[i3=yaCi3; 

d[i3=yaCi3; 

} /* for */ 

*y=yaCns~3; 

for (m=1; m<n; m++) 

C 
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for (i=1; i<=n-m; i++) 



> 



ho=xaCi]-x; 
hp=xaCi+m]“X; 
w=cCi+1]-d[i]; 
if ((den=ho-hp)==0.0) 

ExitOnErrorC'Error in routine POLINT"); 
den=w/den; 
dCi]=hp*den; 
cCi]=ho*den; 

> /* for */ 



*y += (*dy=(2*ns < (n-m) ? cCns+1] : dCns — ])); 
> /* for */ 
f ree((char*)d); 
f ree((char*)c); 
returnC 0 ); 

/ i c ' k ' k ’ K ' k i r k’ k ' k 'k’ k ' k ’k'k / 

/'kirkirk END poLint kkkk/ 



/ 



/kkkkk 

* FUNCTION: windowO 

* 

* This function applies a Blackman window to a vector. 



* Arguments: dataC] 

★ 

* Return value: 0 

* 

* Functions called: none 

* 

* Definitions called: ANSI 

* PI 

* 

* Global variables called: none 

* 

* Significant memory allocation: none 

* 



N 



UNIX 



#if defined ( ANSI ) 

int window(float *data, int N); 

#elif defined ( UNIX ) 

windowC data, N ) 

float *data; 

int N; 

Aendif 

i 

int n; 



for (n=0; n<N; n++) 

dataCn+13=dataCn+13*(0.42+0.5*cos(2.0*PI*(f loat)(n-N/2)/(f loat)(N-D) 
+0.08*cos(4.0*PI*(float)(n-N/2)/(float)(N-1))); 

> /* for */ 
return ( 0 ); 



/kkkkk end window *★★★*/ 

/kkkkk 

* FUNCTION: realftO 

* 

* This function calculates FFT's 

* 

* Arguments: 

* 



dataC3 

isign 



* Return value: 

* 
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fourl ( ) 



★ Functions called: 

★ 

* Definitions called: ANSI 

★ 

* Global variables called: none 

★ 

* Significant memory allocation: none 

* 

irkicirk/ 

Pii defined ( ANSI ) 

int realftCfloat *data, int n, int isign) 

#elif defined ( UNIX ) 
realftCdata, n, isign) 
float *data; 
int n, isign; 

^endif 

int i, i1, 12, i3, i4, n2p3; 

float c1=0.5, c2, hir, h1i, h2r, h2i; 

double wr, wi, wpr, wpi, wtemp, theta; 

theta=3.141592653589793/(double)n; 
if (isign==1) 

c2 = -0.5; 
fourKdata, n, 1); 

> /★ if */ 
else 

C 

c2=0.5; 

theta = -theta; 

> /★ else */ 
wtemp=sin(0.5i*ftheta); 
wpr = -2.0*wtemp*wtemp; 
wpi=sin(theta); 

wr=1 .0+wpr; 

wi=wpi; 

n2p3=2*n+3; 

for (i=2; i<=n/2; i++) 

i 

i4=1+(i3=n2p3-(i2=1+(i1=i+i-1))); 
hi r=c1*(data[i1D+dataCi33); 
h1 i=c1*(data[i2]-data[i43); 
h2r = -c2*(data[i23+dataCi4]); 
h2i=c2*(data[i1 ]-dataCi3]); 
data[i1D=h1 r+wr*h2r-wi*h2i; 
dataCi2]=h1 i+wr*h2i+wi*h2r; 
data[i3]=h1 r-wr*h2r+wi*h2i ; 
data[i4] = -h1i+wr*h2i+wi*h2r; 
wr=(wtemp=wr)*wpr-wi*wpi+wr; 
wi=wi*wpr+wtemp*wpi+wi; 

> /★ for ★/ 
if (isign==1) 

C 

dataCI ]=(h1 r=data[1 ])+dataC23; 
data[2]=h1 r-dataC2]; 

> /* if ★/ 
else 

dataC1]=c1*((h1r=data[1])+dataC2]); 

data[2]=c1*(h1r-dataC2]); 

fourl (data,n,-1); 

> /* else ★/ 
return ( 0 ); 

/tkk'k'k end realft kkirkk/ 
/kkkkkirkirkirkkkkk-k’kk'kirk’k / 

/•kirkirk 



UNIX 
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* FUNCTION; fourlO 

★ 

* This function calculates FFT's 



* 

* Arguments: dataC] 

* isign 

* 

* Return value: 0 

* 

* Functions called: none 

★ 

* Definitions called: ANSI 

* SWAP 

* 

* Global variables called: none 

* 

* Significant memory allocation: none 

* 

tfif defined ( ANSI ) 



int fourl (float *data, int nn, int isign) 

#elif defined ( UNIX ) 

fourKdata, nn, isign) 

float *data; 

int nn, isign; 

#endif 

int n, mmax, m, j, istep, i; 

double wtemp, wr, wpr, wpi, wi, theta; 

float tempr, tempi; 



nn 



UNIX 



n=nn « 1 ; 
i=1; 

for (i=1; i<n; i+=2) 
if (j > 1) 

SWAPCdataCj ],dataCi]); 

SWAP(dataC j+1 ],dataCi+l]); 

> /* for */ 
m=n » 1; 

whi le (m >= 2 && j > m) 

C 

j -= m; 
m »=1; 

> /* while */ 

j += 

y /* for */ 
mmax=2; 

whi le (n > mmax) 

istep=2*mmax; 

theta=6. 2831 853071 7959/ (i si gn*mmax); 

wtemp=sin(0.5*theta); 

wpr = -2.0*wtemp*wtemp; 

wpi=sin(theta); 

wr=1 . 0; 

wi=0.0; 

for (m=1; m<mmax; m+=2) 
i 

for (i=m; i<=n; i+= istep) 
j=i+mmax; 

tempr=wr*dataCj]-wi*dataC j+1]; 
tempi =wr*data[ j+1 ]+wi*data[j ]; 
dataCj]=dataCi]-tempr; 
dataCj+1]=data[i+1]-tempi; 
dataCi] += tempr; 
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data[i+1] += tempi; 

> /★ for ★/ 

wr=(wtemp=wr)*wpr-wi*wpi+wr; 

wi=wi*wpr+wtemp*wpi+wi; 

> /* for */ 
mmax=istep; 

> /★ while ★/ 
return ( 0 ); 

/irkirirk end fourl kkkirk / 

/kkkkk 

★ FUNCTION: vectorO 

* 

This function allocates memory for UNIT OFFSET vectors. 



Arguments: 

Return value: 
Functions called: 
Definitions called: 



length 

*v 

ExitOnErrorO 

ANSI 

none 

v[] 



UNIX 



★ 

★ 

★ 

★ 

★ 

★ 

★ 

★ 

★ 

★ 

★ Global variables called: 

* 

* Significant memory allocation: 

★ 

icirkkk / 

//if defined ( ANSI ) 
float *vector(int length) 

//elif defined ( UNIX ) 
vectorC length) 
int length; 

//endif 

float *v; 



i f ( (v=(f loat*)mal loc( ( length+1 )*si zeof (f loat ) ) )==NULL) 
ExitOnErrorC'Memory allocation failure in vectorO."); 
//if defined ( ANSI ) 
return v; 

//elif defined ( UNIX ) 
return (long int)v; 

//endif 

> /kirkirkirkirkirkirkirkirkkkkkk/ 

/irkirirk END vector irkkkk / 

/irkkirkirkkkkkkkirkkkkkkkk / 



* FUNCTION: matrixO 

* 

★ This function allocates memory for UNIT OFFSET 2-D arrays. 

* 



★ Arguments: row 

* 

* Return value: **m 

* 

* Functions called: ExitOnErrorO 

* 

* Definitions called: ANSI 

★ 

* Global variables called: none 

★ 

* Significant memory allocation: m[][] 

* 

//if defined ( ANSI ) 

float **matrix(int row, int col) 



col 



UNIX 
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#elif defined ( UNIX ) 

matrix(row,col) 

int row, col; 

#endif 

int i; 
float **m; 



if ((m=(f loat**)fna Hoc ((unsigned ) (row+1 )*sizeof (float*) ))==NULL) 
ExitOnErrorC'Allocation failure 1 in matrixO"); 
for (i=1; i<=row; i++) 

if ((m[i]=(f loat*)malloc( (unsigned) (co 1+1 )*sizeof (f loat)) )==NULL) 
ExitOnErrorC'Allocation failure 2 in matrixO"); 

> /* for */ 

#if defined ( ANSI ) 
return m; 

Pelii defined ( UNIX ) 



return (long int)m; 

#endif 

/tirkirk END matrix *****/ 

/ 1 ***** k k k kirk’k’k'kirk'k'k'k'k'k / 

/***** 

* FUNCTION: free^matrixO 

* 

This function deallocates memory from UNIT OFFSET 2-D arrays. 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

★*★*★/ 

#if defined ( ANSI ) 

void f ree_matrix(f loat **m, int row) 

/^elif defined ( UNIX ) 

f ree_matrix(m, row) 

float **m; 

int row; 

#endif 

i 

int i; 



Arguments: m[][] 

Return value: 0 

Functions called: none 

Definitions called: ANSI 

Global variables called: none 

Significant memory allocation: none 



UNIX 



for(i=row; i>=1; i — ) 
f ree( (char*)mCi]); 
f ree((char*)m); 
return ( 0 ); 



/*★★** end free matrix *****/ 



/***** 

* FUNCTION: ExitOnErrorO 

* 

* This function performs an abnormal process termination. 

* 

* Arguments: error_txt[] 

* 

* Return value: none 

* 

* Functions called: none 
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ANSI 



UNIX 



★ 

* Definitions called: 

★ 

* Global variables called: none 

★ 

* Significant memory allocation: none 

★ 



irkirirk/ 

#if defined ( ANSI ) 

void ExitOnErrorCchar error_txtC]) 

#elif defined ( UNIX ) 

ExitOnError (error_txt) 
char error_txtCD; 

#endif 



> 



fprintfCstderr, "Program run-time error ...\n"); 
f printf (stderr, "%s\n",error_txt); 
fprintf (stderr/‘. . .now exiting to system. .. \n"); 
exit(O); 



/krkkirk END ExitOnError kkkirk/ 



B. DECOMMA.C 

This program removes commas and colons from the output of 
the upper tilt sensor of the Heard Island Array. 



* PROGRAM DECOMMA.C vsn 1.0 

* WRITTEN BY: Steven Crocker 

* LAST UPDATE: August 3, 1991 

* 

* This program takes any ASCII file and copies it to a user 

★ defined file. It removes the commas and colons. 
kickicic / 

^define VERSION ANSI 
#define C 25 
^define COMMA 44 
^define COLON 58 
^define SPACE 32 
#include<stdio.h> 

# include <malloc.h> 

#if (VERSION==ANSI) 

#include<stdlib.h> 
void ExitOnErrorCchar error_txtC]); 

#endif 

mainO 

FILE *sacm_fp, *out_fp; 
char c, sacm_f i leCC], out_fileCC3; 

printf ("\n\n\nEnter file name for input SACM data.Xn"); 
scanf ("%s",sacm_f i le); 
printf ("\n\n"); 

printf ("Enter file name for output SACM data.Xn"); 
scanf ("%s"^out_fi le); 

if ((sacm_f p=fopen(sacm_f i le, "rt")> != NULL) i 
if ((out_fp=fopen(out_f i le, "wt")) != NULL) { 
whi le((c=fgetc( sacm_fp )) != EOF) -C 



/★ either ANSI or UNIX */ 
/* max length for file names k/ 
/★ ASCII codes */ 
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if((c==COMMA) 11 (c==COLON)) c=SPACE; 
fputcCc, out_fp); 

> /* whi Le */ 

> /* if ★/ 

else ExitOnErrorC'Error opening output data file"); 

> /★ if */ 

else ExitOnErrorC'Error opening input data file"); 
f close(sacni_fp); 
fclose(out_fp); 
exit(O); 

> 

#if (VERSI0N==ANSI) 

void ExitOnError (char error_txtCD) 

#elif (VERSI0N==UN1X) 

ExitOnError (error_txt) 
char error_txtC]; 

#endif 

i 

fprintfCstderr, "Program run-time error ...\n"); 
fprintf(stderr,"%s\n",error_txt); 
fprintfCstderr, .now exiting to system. .. \n"); 
exit n ); 

> 



C . SACMl . C 



This program reads the data output by the upper tilt 
sensor and formats it for use in the beamformer. It 



calculates 60 second averages for all data fields, converts 
pressure to depth and calculates tilt direction based on 
current velocity. 



/AAAAA 

* PROGRAM SACH1.C vsn 1.0 

* WRITTEN BY; Steven Crocker 

* UST UPDATE: August 6, 1991 
A 

* This program takes SACM data from the Heard Island West Coast 

* Array (upper instrument package) and condenses it. 

* 

* Pressure is converted to depth. 

■k 

* The conductivity is not processed. All values output are 60 

* second averages of the input. 

* 



irkirkk/ 

#def ine 


VERSION ANSI 


/* 


either ANSI or UNIX 


*/ 


//define 


SELECT OUTPUT OFF 


/* 


ON or OFF 


*/ 


//define 


C 25 


/* 


max length for file names 


*/ 


//define 


inBUFFER 1200 


/* 


input buffer size 


*/ 


//define 


outBUFFER 50 


/* output buffer size 


★/ 


//define 

//define 


PI 3.14159265359 
RADIAN 57.2957795131 








#include<stdio. h> 
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^include <malloc.h> 
#include<math. h> 



#if (VERSION==ANSI) 

#include<stdLib.h> 

void ExitOnErrorCchar error_txt[]); 

#endif 



mainO 



FILE *sacm_fp, *out_fp; 

char c, sacm_f i leCC], out_file[C]; 

int i, j, count, loop_count, inBUFFER^fuLl, out_flag[6]; 

float t1, avgTilt, avgPress, avgVN, avgVE, avgTemp, avgAngle, 
avgVelocity, junk; 

static float inTi ItCinBUFFER], inPressCinBUFFER], inVNCinBUFFER], 
inVECinBUFFER], inTiineCinBUFFER], inTemplinBUFFER], 
outTimeCoutBUFFER], outTi It CoutBUFFER], 
out PressEout BUFFER], outAng LeCout BUFFER], 
outVelocity [out BUFFER], out TempCout BUFFER]; 

printf ("\n\n\nEnter file name for SACM data.Xn"); 
scanf ( "%s" , sacm_f i le) ; 
printf ("\n\n"); 

printf ("Enter file name for output. \n"); 
scanf ("%s",out_fi Le); 

printf ( "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n" ); 
printf ("The following pertains to output file format . \n\n"); 

printf ("Units of measure: Time index is minutes. \n"); 

printfC Angles in degreesXn"); 

printfC Compass direction (TRUE)\n"); 

printfC Velocity in meters/secondXn"); 

printfC Temperature in degrees C\n"); 

printfC Depth in meters\n\n"); 

printf ("NOTE: Output column order is same as above. \n\n\n"); 

for ( i=0; i<6; i++) out_f lag[i]=1; 

#if (SELECT_OUTPUT==ON) 

printf ("Include TIME? (y or n)\n"); 

scanf("%s",&c); 

if(c==89 II c==121) out_flag[0]=1; 
else out_f lag[0]=0; 
printfC\n\n"); 

printf ("Include TILT? (y or n)\n"); 
scanf ("%s",&c); 

if(c==89 !i c==121) out_f lag[1]=1; 
else out_f lagC1]=0; 
printf ("\n\n"); 

printfClnclude CURRENT DIRECTION? (y or n)\n"); 
scanf C%s",&c); 

if(c==89 11 c==121) out_f lag[2]=1; 
else out_f Lag[2]=0; 
printfC\n\n"); 

printfClnclude CURRENT VELOCITY? (y or n)\n"); 
scanfC%s",&c); 

if(c==89 i; c==121) out_flag[3]=l; 
else out_f lagC3]=0; 
printfC \n\n"); 

printfClnclude TEMPERATURE? (y or n)\n"); 
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scant ("%s",&c); 

if(c==89 11 c==121) out_flagC4]=1; 
else out_f LagC4]=0; 
printf ("\n\n"); 

printf ("Include DEPTH? (y or n)\n"); 
scant ("%s",&c); 

it(c==89 11 c==121) out_t lagC5]=1; 
else out_t lagC5]=0; 

#endi t 

inBUFFER_tuU=inBUFFER; 

loop_count=0; 

it ((sacm_tp=topen(sacm_ti Le, "rt")) != NULL) ; 
else ExitOnErrorC'Error opening input data tile"); 
it ( (out_tp=topen(out_ti le, "wt")) != NULL) ; 
else ExitOnErrorC'Error opening output data tile"); 

while (1) 

tor (i=0; i<inBUFFER; i++) 

it (tscant(sacm_tp, "%t %t %t %t %t %t %t\n", 

&inTimeCi],&inVNCi3, SinVECiD, SinTempCi], 
SinTiltCi]/ SinPressCi], Sjunk) != EOF) ; 
else inBUFFER_tull=i-2; 

} /* tor */ 
j=0; 

tor (i=0; i<inBUFFER tull;i++) 

i 

t1=inTimeCi3; 

avgTi lt=0.0; 

avgPress=0.0; 

avgVN=0.0; 

avgVE=0.0; 

avgTemp=0.0; 

count=0; 

whi le((inTimeCi]<t1+60.0) && (i<inBUFFER tulD) 

i 

avgTi lt=avgTi It + inTUtCi]/10.0; 
avgPress=avgPress + inPressLi]*1000.0; 
avgVN=avgVN + inVNCi 3/1000. 0; 
avgVE=avgVE + inVE[i3/1000.0; 
avgTemp=avgTemp + inTempCi3/100.0; 
i++; 

count++; 

> /* whi le */ 

avgTi lt=avgTi lt/(tloat) count; 
avgPress=avgPress/(t loat)count ; 
avgTemp=avgTemp/(t loat)count; 
avgVN=avgVN/(t loat) count; 
avgVE=avgVE / ( t loat ) count ; 
avgVelocity=sqrt (avgVN*avgVN+avgVE*avgVE); 
avgAngle=atan2(avgVE, avgVN)*RADIAN; 
it (avgAngle<=0.0) avgAngle=360.0+avgAngle; 
outTimeC j3=(t loat) ( j+1)+outBUFFER*loop_count; 
outTi ItCj 3=avgTi It; 

out Pr essC j 3= (avgPress-1 01 325. 0)7(9.80665*1 026.0); 
outTempC j3=avgTemp; 
outVe loci ty[j3=avgVe loci ty; 
outAngleC j3=avgAngle; 

j++; 

} /* tor */ 



it (inBUFFER_tull != inBUFFER) j— ; 
tor (i=0; i<j; i++) 

C 

it (out_t lag[03) tprintt (out_tp, "%6.0t", outTimeCi3) 
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if (out__f LagCI]) fprintf (out_fp, 
if (out~f LagC23) fprintf (out^fp, 
if (out_f Lag[33) fprintf (out_f p, 
if (out_f LagC43) fprintf (out_fp, 
if (out_f LagC53) f printf (out_fp, 
fprintf (out_fp, "\n"); 

> /* for */ 



%8.5f, outTiltCi3); 
%8.3f", outAngLeCi3); 
%8.6f", outVeLocityCi3); 
%8.5f", outTempCi3); 
%8.0f“, outPressCi3); 



if (inBUFFER^fuLL != inBUFFER) exit(O); 
loop_count++; 

> /* while */ 
exit(l); 



#if (VERSI0N==ANSI) 

void ExitOnErrorCchar error__txtC3) 

#elif (VERSION==UNIX) 

ExitOnError(error_txt) 
char error_txtC3; 

#endif 

fprintf (stderr, •'Program run-time error ...\n"); 
fprintf (stderr, "%s\n",error_txt); 
fprintf (stderr/'. . .now exiting to system...\n"); 
ex i t ( 1) ; 

> 



D . S ACM2 . C 



This program takes the output of sacml.c as input. It 
locates and copies a user defined subset of the tilt data for 



use in beamforming. 



/irk-kirk 

* PROGRAM SACM2.C vsn 1.0 

* WRITTEN BY: Steven Crocker 

* LAST UPDATE: August 6, 1991 

* 

* This program takes data processed by SACM1.C and cuts 

* a user defined segment from it. The segment is retained 

* as the output data file. The input file is not affected. 

* The output time base may either be normalized to 1 or may 

* retain the original values. 

k 

* The option is given to accept a default output format. This 

* output format coincides with the required input format to the 

* time domain modal beamformer. 

* 

* The output elements are selectable if desired. 

^define VERSION ANSI 
^define C 25 
#include<stdio. h> 

^include <malloc.h> 

#if (VERSION==ANSI) 

#include<stdlib.h> 
void ExitOnErrorCchar error_txtC3); 

^endif 



/* either ANSI or UNIX */ 

/* max length for file names */ 
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ma i n ( ) 
C 



FILE *input_fp, *out_fp; 

char c. Default, input_f i le[C], out_file[C]; 

int i=0, out1, outN, out_flag[7]; 

float Tilt, Depth, Time, Temp, Angle, Velocity; 

printf ("\n\n\nEnter file name for input data.Xn"); 
scanf (“%s", input_f i le) ; 

if ( ( input_fp=fopen( input_f i le, "rt")) != NULL) ; 
else ExitOnErrorC "Error opening input file"); 
printf("\n\n"); 

printf ("Enter file name for output. \n"); 
scanf ("%s",out_f i le); 

if ( (out_fp=fopen(out_f i le, "wt")) != NULL) ; 
else ExitOnErrorC'Error opening output file"); 

printf ("\n\n\nDo you want the default output conf iguration?\n"); 
printf ("This option formats the output for use in the beamformer . \n"); 
scanf("%s", Default); 

printf ( "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n" ); 
if(Default != ‘n’ &S Default != *N') 



printf ("The following pertains to output file format , \n\n"); 
printf ("Units of measure: Time Index is minutes. \n"); 



printf ("NOTE: Output column order is same as above. \n\n\n"); 

printf ("Include TIME? (y or n)\n"); 
scanf ("%s",c); 

if(c==89 II c==121) out_flag[0]=1; 
else out_f lag[0D=0; 
printf ("\n\n"); 

printf ("Include TILT? (y or n)\n"); 
scanf C7»s",c); 

if(c==89 II c==121) out_f lag[1]=1; 
else out_f lag[1]=0; 
printf ("\n\n"); 

printfC'Include CURRENT DIRECTION? (y or n)\n"); 
scanf ("%s",c); 

if(c==89 II c==121) out_flag[2]=1; 
else out_f lag[2]=0; 
printf ("\n\n"); 

printfC'Include CURRENT SPEED? (y or n)\n"); 
scanf ("%s",c); 

if(c==89 II c==121) out_flag[3]=1; 
else out_f lagC3]=0; 
printf ("\n\n"); 

printfC'Include TEMPERATURE? (y or n)\n"); 
scanf ("%s",c); 

if(c==89 II c==121) out_flag[4]=1; 
else out_f lag[4]=0; 
printf ("\n\n"); 

printfC'Include DEPTH? (y or n)\n"); 
scanf ("%s",c); 

if(c==89 11 c==121) out_flag[5]=1; 
else out_f lagC53=0; 



printf (" 
printf (" 
printf (" 
printf (" 
printf (" 



Angles in degreesXn"); 



Compass direction (TRUE)\n"); 
Velocity in meters/secondXn"); 
Temperature in degrees CXn"); 
Depth in metersXnXn"); 



> /* if */ 
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else 

out_f lagC0]=0; 
out_f LagC1]=1; 
out_f Lag[2]=1; 
out_f lag[3]=0; 
out_f lag[4]=0; 
out flag[5]=1; 

> /* else */ 

printf ("\n\nWhat time index of the input file should be the firstXn"); 

printf ("element of the output file?\n"); 

scanf("%i",&out1); 

printf("\n\n"); 

printf ("How many elements should the output file contain?\n"); 
scanf ("%i",&outN); 
p^intf ("\n\n"); 

if (Default != 'n' && Default != 'N') 

printf ("Should the time base be normalized to begin at t=1?\n"); 
scanf ("%s", c); 

if(c==89 ;i c==121) out_f lag[6]=1; 
else out_f lag[6]=0; 
printf ("\n\n\n"); 

> /* if */ 



> 



while (1) 
i 



if (fscanf (input_fp, "%f %f %f %f %f %f \n", 

STime, &Tilt, SAngle, ^Velocity, STemp, SDepth) != EOF) 

£ 

if((Time >= (float)outl) && (Time < (f loat) (out1+outN)) ) 



if(out_f lagCO] && lout^f lag[6]) 
fprintf (out_fp/'%f ", Time); 

else if (out_flagC0] && out_f lag[6!l) 

fprintf(out_fp, "%f", (f loat)(i+D); 
if (out_f lag[1 ]) fprintf (out_fp, 
if (out_f lag[2]) fprintf (out_fp, 
if (out_f lag[3]) fprintf (out_fp, 
if (out_f lag[4]) fprintf (out_fp, 
if (out_f lag[5]) fprintf (out_fp, 
fprintf (out_fp, "\n"); 
i++ ; 

> /* if */ 

else if(Time >= (f loat) (out1+outN)) 



y.i", 


Tilt); 


Xf", 


Angle); 


Xf, 


Velocity); 


7.i\ 


Temp); 


Xf", 


Depth); 



f close( input_f p); 
f close(out_fp); 
exit(O); 

> /* else if */ 

> /* if */ 

else 

i 



f close(input_fp); 
f close(out_fp); 

Exi tOnErrorC'EOF encountered in input data file"); 
> /* else */ 

> /* while */ 
ex i t ( 1 ) ; 



#if (VERSI0N==ANS1) 

void Exi tOnError(char error_txt[D) 

#elif (VERSION==UNIX) 

Exi tOnError(error txt) 
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char error^txtCD; 

#endif 

fprintfCstderr, "Program run-time error ...\n"); 
fprintf (stderr, "%s\n",error_txt ); 
f printf (stderr, " . . . now exiting to system. .. \n"); 
exit(l); 

> 



E. ARRAYTEST.C 



This program performs various statistical tests on each 
element of the acoustic array. It is useful for isolating 
hydrophones which have failed during the experiment. The 
statistical tests are an adaptation of those found in 



Reference 11. 



/irkirk 



* PROGRAM; ARRAYTST vsn 1.0 

* WRITTEN BY: Steven Crocker 

* LAST UPDATE: October 21, 1991 

* 








★ 

irkirk/ 


^define UNIX VERSION 


/* 


either ANSI or UNIX 


*/ 


/(/define CHANNELS 32 


/★ 


number of hydrophones on array 


*/ 


/(/define F SAMPLE 228 


/* 


sampling frequency 


*/ 


//define BUFFER TIME 60 


/* 


duration of time averages 


*/ 



#include<stdio.h> 

U i nc Lude<ma I loc . h> 
#incLude<math.h> 



#if defined ( ANSI ) 

#incLude<f Loat .h> 

#incLude<stdlib.h> 
int getlnput(void); 

int momentsCf Loat *data, int n, float *ave, float *ave2, float *adev, 
float *sdev, float *svar, float *skew, float *curt) 
float *vector(int Length); 
float **matrix(int row, int col); 
int f ree_matrix(f loat **m, int row); 
void ExitOnErrorCchar error_txt[]); 

#elif defined ( UNIX ) 
int vectorO, 
matrixO, 
getInputO, 

momentsO, 
f ree_matrix(), 

ExitOnErrorO; 

#endif 

/* Global Variable Declarations */ 
float **inSOUND; 
int f irstBuffer=1; 

mainO 

int i, 3=1, k, flagCS]; 
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float ave, ave2, adev, sdev, svar, skew, curt, *temp; 
char f i leNameCSO], key; 

FILE *fp1, *fp2, *fp3, *fp4, *fp5, *fp6, *fp7; 

inSOUND=(f Loat **)matrix(CHANNELS, BUFFER_TIME*F_SAMPLE); 
temp=(f Loat *)vector(BUFFER_TIME*F_$AMPLE); 

for (i=0;i<=7;i++) flagCi3=0; 

pn’ntf ("Answer y or n to the following output options. \n"); 
printf ("Average value ( over %i seconds ): ", BUFFER_TIME); 
$canf("%s", ftkey); 
if(key==89 H key==121) 
i 

flag[13=1; 

printf ("Enter file name: "); 
scanf("%s", filename); 
printf ("\n\n"); 

if ((fp1=fopen(f i leName, "wt")) == NULL) 

ExitOnErrorC'Error opening file."); 

> /* if ★/ 

printf ("Mean squared value ( over %i seconds ): ", BUFFER_TIME) 
scanf("%s", &key); 
if(key==89 11 key==121) 

< 

flag[23=1; 

printf ("Enter file name: "); 
scanf("%s", filename); 
printf ("\n\n"); 

if ( (fp2=fopen(fi leName, "wt")) == NULL) 

ExitOnErrorC'Error opening file."); 

> /* if */ 

printf ("Average deviation ( over %i seconds ): ", BUFFER_TIME); 
scanf("%s", gkey); 
if(key==89 H key==121) 
i 

flagC33=1; 

printf ("Enter file name: "); 
scanf("%s", filename); 
printf ("\n\n"); 

if ( (fp3=fopen(fi leName, "wt")) == NULL) 

ExitOnErrorC'Error opening file."); 

> /* if */ 

printf ("Standard deviation ( over %i seconds ): ", BUFFER_TIME) 
scanf("%s", gkey); 
if(key==89 I I key==121) 

C 

flag[43=1; 

printf ("Enter file name: "); 
scanf("%s", filename); 
printf ("\n\n"); 

if ((fp4=fopen(fi leName, "wt")) == NULL) 

ExitOnErrorC'Error opening file."); 

> /* if */ 

printf ("Variance ( over %i seconds ): ", BUFFER_TIME); 
scanf("%s", gkey); 
if(key==89 H key==121) 

flag[53=1; 

printf ("Enter file name: "); 
scanf("%s", filename); 
printf ("\n\n"); 

if ((fp5=fopen(fi leName, "wt")) == NULL) 

ExitOnErrorC'Error opening file."); 

> /* if ★/ 
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BUFFER_TIME); 



printf ("Skewness ( over %i seconds ): ", 
scanf("%s", &key); 
if(key==89 H key==121) 

C 

f LagC6]=1; 

printf ("Enter file name: "); 
scanf("%s", fileName); 
printf ("\n\n"); 
if ((fp6=fopen(f i LeName, "wt")) == NULL) 
ExitOnErrorC'Error opening file."); 

> /* if */ 



printf ("Kurtosis ( over %i seconds ): ", BUFFER__TIME); 
scanf("%s", &key); 
if(key==89 jj key==121) 
i 

flag[7]=1; 

printf ("Enter file name: "); 
scanf("%s", fileName); 
printf ("\n\n"); 

if ((fp7=fopen(fileName, "wt")) == NULL) 
ExitOnErrorC'Error opening file."); 

> /* if */ 

while (j<120) 

get Input (); 

if (flagCI]) fprintf (fp1,"%i", j ); 
if (flagC2]) fprintf(fp2,"%i",j); 
if (flag[3]) fprintf(fp3,"%i",j); 
if (flag[4]) fprintf (fp4, "%i ",j ); 
if (flag[5]) fprintf (fp5,"%i",j); 
if (flagC6]) fprintf (fp6,"%i",j); 
if (flag[7]) fprintf(fp7,'7ci",j); 

for (i=1;i<=CHANNELS;i++) 
i 

for(k=1;k<=BUFFER_TIME*F_SAMPLE;k++) 

temp[k]=inSOUND[i][k]; 

moments(temp,BUFFER_TIME*F_SAMPLE,&ave,&ave2, 

&adev,&sdev,&svar,&skew,&curt); 



if 


(flagCI]) fprintf(fp1. 


" %f 


", ave); 


if 


(flag[2]) fprintf(fp2. 


" %f 


", ave2); 


if 


(flag[3]) fprintf(fp3. 


" %f 


", adev); 


if 


(flagCA]) fprintf(fp4. 


" %f 


", sdev); 


if 


(flagCS]) fprintf(fp5. 


" 


", svar); 


if 


(flagC6]) fprintf(fp6. 


" %f 


", skew); 


if 


(flagC7]) fprintf(fp7. 


" %f 


", curt); 



> /* for */ 

if (flagCI]) fprintf(fp1,"\n"); 
if (flagC2]) fprintf (fp2,"\n"); 
if (flagC3]) fprintf (fp3,"\n"); 
if (flagCA]) fprintf (fp4,"\n"); 
if (flagCS]) fprintf (fp5,"\n"); 
if (flagC63) fprintf (fp6, "\n"); 
if (flag[7]) fprintf(fp7/'\n"); 
f irstBuffer=0; 

j++; 

> /* whi le */ 



exit(O); 




/irktirk 

* FUNCTION: getInputO 
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* 

* 

* 

* 

* 

* 

* 

★ 

* 

* 

★ 

★ 

★ 

* 

★ 

★ 

* 

* 



This function handles all acoustic input. It also provides one of 
two normal process terminations available in the program. 



Arguments: 


none 




Return value: 


0 




Functions called: 


vectorO 


ExitOnErrorO 


Definitions called: 


ANSI 


UNIX 




F^SAMPLE 

BUFFER_TIME 


CHANNELS 


Global variables called: 


inSOUNDCDC] 


fi rstBuf fer 


Significant memory allocation: 


diskBuf ferC] 





* * ** A / 

#if defined ( ANSI ) 
int getinput (void) 

#elif defined ( UNIX ) 
getInputO 
# end if 

int i, j, items, buffer; 
float ★diskBuffer; 
char f i leNameCSO]; 
static FILE *fpStatic; 



if (f i rstBuf f er) 

printf ("Enter file name for input acoustic data: "); 



if ((scanf ("%s", f i leName) )==EOF) 

Exi tOnErrorC'Fatal error in scanfO"); 



printf("\n\n"); 

if ( (fpStatic=fopen(fi leName, "rb")) == NULL) 

ExitOnErrorC'Error opening INPUT ACOUSTIC data file."); 
> /* if ★/ 



/★ Memory Allocation */ 

buffer=BUFFER_TIME*F_SAMPLE*CHANNELS; 
diskBuffer=(f loat*) vector (buffer); 



items=fread((char*)(diskBuffer+1), sizeof (float), buffer, fpStatic); 
if (items==buf fer) /*continue*/; 
else if (f error(fpStatic) != 0) 

ExitOnErrorC'Error encountered while reading input acoustic data") 
else if (feof (fpStatic) != 0) 



printf (" 
printf ("\t 
printf (" 
fclose(fpStatic); 

exi t (0); 
/★ else if */ 



End of File reached: EXECUTION COMPLETEXn"); 



else ExitOnErrorC'Unknown error handling acoustic input file."); 



for (i=1; i<=BUFFER TIME*F_SAMPLE; i++) 

for (j=1; j<=CHANNELS; j++) 

inSOUNDCjDCi] = disk3ufferCCHANNELS*(i-1)+j]; 

> /* for ★/ 



/* Deallocate Memory ★/ 
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f ree((char*)diskBuf fer); 
return ( 0 ); 

/*★*** end get Input 

/ ★★ k - kk -ki r kic- k ' k / 

/tirkirk 



* 

■jf 


FUNCTION: momentsO 






* 


This function calculates the 


mean, mean squared value, average 


* 


deviation, standard deviation 


, variance, skewness and kurtosis 


k 


of a data vector. 






* 


Arguments: 


data 


n 


* 




ave 


ave2 


★ 




adev 


sdev 


★ 




svar 


skew 


★ 




curt 




★ 








* 


Return value: 


0 




★ 


Functions called: 


ExitOnErrorO 




* 


Definitions called: 


ANSI 


UNIX 


★ 


Global variables called: 


none 




* 

* 


Significant memory allocation 


none 




#if defined ( ANSI ) 

int moments(f loat ★data, int n. 


float ★ave, float 


★ave2, float ★adev, 



float *sdev, float *svar, float *skew, float *curt) 
#elif defined ( UNIX ) 

moments (data, n,ave,ave2,adev,sdev,svar, skew, curt) 
int n; 

float *data, *ave, *ave2, *adev, *sdev, *svar, *skew, *curt; 
#endif 

int j; 
float s,p; 



if (n<=1) ExitOnErrorC'n must be at least 2 in momentO*'); 
s=0.0; 

*ave2=0.0; 

for ( j=1; j<=n; j++) 

s += data[j]; 

*ave2 += data[j]*data[j ]; 

> /* for */ 

*ave=s/n; 

*ave2 /= n; 

*adev=(*svar)=(*skew)=(*curt)=0.0; 
for ( j=1; j<=n; j++) 

★adev += fabs(s=dataC j]“(*a\/e) ); 

★svar += (p=s*s); 

★skew += (p ★= s); 

★curt += (p ★= s); 

} /* for ★/ 

★adev /= n; 

★svar /= (n-1); 

★sdev=sqrt (★svar) ; 
if (★svar) 

< 

★skew /= (n^(^svar)^(^sdev) ); 
★curt=(★curt)/(n★(★svar)★(*svar))“3.0; 

> /★ if */ 

else ExitOnErrorC'No skew/kurtosis when variance = 0 (in momentO ) "); 
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returnC 0 ); 

t «-i i I t * I 1 1 I i 1 / 

/ RXxwTKTTTirTirxiciritxxxx / 

/★★★★* end moment irtrkirk/ 

/■kirkirk 

* FUNCTION: vectorO 

★ 

* This function allocates memory for UNIT OFFSET vectors. 

★ 

* Arguments: length 

★ 

* Return value: *v 

★ 

* Functions called: ExitOnErrorO 

* 

★ Definitions called: ANSI UNIX 

★ 

* Global variables called: none 

★ 

* Significant memory allocation: vC] 

★ 



#if defined ( ANSI ) 
float *vector(int length) 
#elif defined ( UNIX ) 
vectorC length) 
int length; 

#endif 

i 

float *v; 



if ((v=(f loat*)malloc(( length+1)*sizeof (f loat) ))==NULL) 
ExitOnError(‘'Memory allocation failure in vectorO."); 
#if defined ( ANSI ) 
return v; 

#elif defined ( UNIX ) 



return (long int)v; 

#endif 

> /irkkirkkickkk’k’kkkkkirir’kirkk / 

/irkkirk END vector kkkkk / 

/kkkkkkkkkkkkkkkkkkkkkk / 

/★★★★★ 

* FUNCTION: matrixO 

* 

* This function allocates memory for UNIT OFFSET 2-D arrays. 

★ 



* Arguments: row 

★ 

* Return value: **m 

* 

* Functions called: ExitOnErrorO 

k 

* Definitions called: ANSI 

* 

* Global variables called: none 

★ 

* Significant memory allocation: m[]E] 

★ 

kkkkk f 

#if defined ( ANSI ) 

float *^matrix(int row, int col) 

#elif defined ( UNIX ) 

matrix(row,col) 

int row, col; 

#endif 



col 



UNIX 



int i; 
float **m; 
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if ( (m=( float**) mat Loc( (unsigned) (row+1 )*sizeof (f loat*)))==NULL) 
ExitOnErrorC'Allocation failure 1 in tnatrixO"); 
for (i=1; i<=row; i++) 

if ( (m[i]= (float* )ma I loc( (unsigned) ( CO 1+1 )*si zeof (f loat )) )==NULL) 
ExitOnErrorC'Allocation failure 2 in matrixO"); 

} /* for */ 

#if defined ( ANSI ) 
return m; 

#elif defined ( UNIX ) 
return (long int)m; 

#endif 

> /*j^»*j>r * * A*»**A* * * ******/ 

/ir kk ** end matrix **★**/ 

/****» * * *****» ** *** *★* * / 

/irklrirk 

* FUNCTION: f ree^matrixO 

★ 

* This function deallocates memory from UNIT OFFSET 2-D arrays. 

★ 

* Arguments: m[][] row 

★ 



k 

k 


Return value: 


0 


k 

k 


Functions called: 


none 


k 

k 


Definitions called: 


ANSI 


k 

k 


Global variables called: 


none 


k 

★ 


Significant memory allocation: 


none 



icirkirk/ 

#if defined ( ANSI ) 

void f ree_matrix(f loat **m, int row) 

#elif defined ( UNIX ) 

f ree_matrix(m, row) 

float **m; 

int row; 

#endif 

C 

int i; 



for(i=row; i>=1; i — ) 
f ^ee((char*)m[i]); 
f ree((char*)m); 
return( 0 ); 



/kirktk END free matrix *****/ 



/ kk kkk 

* FUNCTION: ExitOnErrorO 

★ 

* This function performs an abnorma 

★ 

* Arguments: 

★ 

* Return value: 

★ 

* Functions called: 
k 

* Definitions called: 

* 

* Global variables called: 

★ 

* Significant memory allocation: 

★ 

■kkk k k/ 



I process termination. 

error_txtC!l 

none 

none 

ANSI UNIX 

none 

none 
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#if defined ( ANSI ) 

void Exi tOnErrorCchar error_txt[]) 

#elif defined ( UNIX ) 

ExitOnError(error_txt) 
char error_txtCD; 

^endif 

C 

fprintf (stderr/'Program run-time error ...\n"); 
fprintf (stderr, "%s\n",error_txt); 
fprintf (stderr,". . .now exiting to system. .. \n"); 
exit(O); 

end ExitOnError 
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